home *** CD-ROM | disk | FTP | other *** search
/ STraTOS 1997 April & May / STraTOS 1 - 1997 April & May.iso / CD01 / INTERNET / SITES / LITTLE / P3SRC.ZIP / ATARI / BLOB.C < prev    next >
Encoding:
C/C++ Source or Header  |  1997-01-18  |  75.6 KB  |  3,806 lines
  1. /****************************************************************************
  2. *                   blob.c
  3. *
  4. *  This module implements functions that manipulate blobs.
  5. *
  6. *  The original file was written by Alexander Enzmann.
  7. *  He wrote the code for blobs and generously provided us these enhancements.
  8. *
  9. *  Modifications and enhancements by Dieter Bayer [DB].
  10. *
  11. *  from Persistence of Vision(tm) Ray Tracer
  12. *  Copyright 1996 Persistence of Vision Team
  13. *---------------------------------------------------------------------------
  14. *  NOTICE: This source code file is provided so that users may experiment
  15. *  with enhancements to POV-Ray and to port the software to platforms other
  16. *  than those supported by the POV-Ray Team.  There are strict rules under
  17. *  which you are permitted to use this file.  The rules are in the file
  18. *  named POVLEGAL.DOC which should be distributed with this file. If
  19. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  20. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  21. *  Forum.  The latest version of POV-Ray may be found there as well.
  22. *
  23. * This program is based on the popular DKB raytracer version 2.12.
  24. * DKBTrace was originally written by David K. Buck.
  25. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  26. *
  27. *****************************************************************************/
  28.  
  29. /****************************************************************************
  30. *
  31. *  Explanation:
  32. *
  33. *    -
  34. *
  35. *  Syntax:
  36. *
  37. *    blob
  38. *    {
  39. *      threshold THRESHOLD_VALUE
  40. *
  41. *      component STRENGTH, RADIUS, <CENTER>
  42. *
  43. *      sphere { <CENTER>, RADIUS, [strength] STRENGTH
  44. *        [ translate <VECTOR> ]
  45. *        [ rotate <VECTOR> ]
  46. *        [ scale <VECTOR> ]
  47. *        [ finish { ... } ]
  48. *        [ pigment { ... } ]
  49. *        [ tnormal { ... } ]
  50. *        [ texture { ... } ]
  51. *      }
  52. *
  53. *      cylinder { <END1>, <END2>, RADIUS, [strength] STRENGTH
  54. *        [ translate <VECTOR> ]
  55. *        [ rotate <VECTOR> ]
  56. *        [ scale <VECTOR> ]
  57. *        [ finish { ... } ]
  58. *        [ pigment { ... } ]
  59. *        [ tnormal { ... } ]
  60. *        [ texture { ... } ]
  61. *      }
  62. *
  63. *      [ sturm ]
  64. *      [ hierarchy FLAG ]
  65. *    }
  66. *
  67. *  ---
  68. *
  69. *  Jul 1994 : Most functions rewritten, bounding hierarchy added. [DB]
  70. *
  71. *  Aug 1994 : Cylindrical blobs added. [DB]
  72. *
  73. *  Sep 1994 : Multi-texturing added (each component can have its own texture).
  74. *             Translation, rotation and scaling of each component added. [DB]
  75. *
  76. *  Oct 1994 : Adopted the method for the bounding slab creation to build the
  77. *             bounding sphere hierarchy of the blob to get a much better
  78. *             hierarchy. Improved bounding sphere calculation for tighter
  79. *             bounds. [DB]
  80. *
  81. *  Dec 1994 : Added code for dynamic blob queue allocation. [DB]
  82. *
  83. *  Feb 1995 : Moved bounding sphere stuff into a seperate file. [DB]
  84. *
  85. *****************************************************************************/
  86.  
  87. #include "frame.h"
  88. #include "povray.h"
  89. #include "vector.h"
  90. #include "povproto.h"
  91. #include "blob.h"
  92. #include "bbox.h"
  93. #include "lighting.h"
  94. #include "matrices.h"
  95. #include "objects.h"
  96. #include "polysolv.h"
  97. #include "texture.h"
  98.  
  99.  
  100.  
  101. /*****************************************************************************
  102. * Local preprocessor defines
  103. ******************************************************************************/
  104.  
  105. /* Minimal intersection depth for a valid intersection. */
  106.  
  107. #define DEPTH_TOLERANCE 1.0e-2
  108.  
  109. /* Tolerance for inside test. */
  110.  
  111. #define INSIDE_TOLERANCE 1.0e-6
  112.  
  113. /* Tolerance used for order reduction during root finding. */
  114.  
  115. #define ROOT_TOLERANCE 1.0e-3
  116.  
  117. /* Ray enters/exits a component. */
  118.  
  119. #define ENTERING 0
  120. #define EXITING  1
  121.  
  122.  
  123.  
  124. /*****************************************************************************
  125. * Local typedefs
  126. ******************************************************************************/
  127.  
  128. typedef struct
  129. {
  130.   int type;
  131.   DBL bound;
  132.   BLOB_ELEMENT *Element;
  133. } Blob_Interval;
  134.  
  135.  
  136.  
  137. /*****************************************************************************
  138. * Static functions
  139. ******************************************************************************/
  140.  
  141. static void element_normal PARAMS((VECTOR Result, VECTOR P, BLOB_ELEMENT *Element));
  142. static int intersect_element PARAMS((VECTOR P, VECTOR D, BLOB_ELEMENT *Element, DBL mindist, DBL *t0, DBL *t1));
  143. static void insert_hit PARAMS((BLOB_ELEMENT *Element, DBL t0, DBL t1, Blob_Interval *intervals, int *cnt));
  144. static int determine_influences PARAMS((VECTOR P, VECTOR D, BLOB *Blob, DBL mindist, Blob_Interval *intervals));
  145. static DBL calculate_field_value PARAMS((BLOB *Blob, VECTOR P));
  146. static DBL calculate_element_field PARAMS((BLOB_ELEMENT *Element, VECTOR P));
  147.  
  148. static int intersect_cylinder PARAMS((BLOB_ELEMENT *Element, VECTOR P, VECTOR D, DBL mindist, DBL *tmin, DBL *tmax));
  149. static int intersect_hemisphere PARAMS((BLOB_ELEMENT *Element, VECTOR P, VECTOR D, DBL mindist, DBL *tmin, DBL *tmax));
  150. static int intersect_sphere PARAMS((BLOB_ELEMENT *Element, VECTOR P, VECTOR D, DBL mindist, DBL *tmin, DBL *tmax));
  151. static int intersect_ellipsoid PARAMS((BLOB_ELEMENT *Element, VECTOR P, VECTOR D, DBL mindist, DBL *tmin, DBL *tmax));
  152.  
  153. static void get_element_bounding_sphere PARAMS((BLOB_ELEMENT *Element, VECTOR Center, DBL *Radius2));
  154. static void build_bounding_hierarchy PARAMS((BLOB *Blob));
  155.  
  156. static void init_blob_element PARAMS((BLOB_ELEMENT *Element));
  157. static void determine_element_texture PARAMS((BLOB *Blob,
  158.   BLOB_ELEMENT *Element, TEXTURE *Texture, VECTOR P, int *Count,
  159.   TEXTURE **Textures, DBL *Weights));
  160.  
  161. static void insert_node PARAMS((BSPHERE_TREE *Node, int *size));
  162.  
  163. static int  All_Blob_Intersections PARAMS((OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack));
  164. static int  Inside_Blob PARAMS((VECTOR point, OBJECT *Object));
  165. static void Blob_Normal PARAMS((VECTOR Result, OBJECT *Object, INTERSECTION *Inter));
  166. static void *Copy_Blob PARAMS((OBJECT *Object));
  167. static void Translate_Blob PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  168. static void Rotate_Blob PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  169. static void Scale_Blob PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  170. static void Invert_Blob PARAMS((OBJECT *Object));
  171. static void Transform_Blob PARAMS((OBJECT *Object, TRANSFORM *Trans));
  172. static void Destroy_Blob PARAMS((OBJECT *Object));
  173. static void Compute_Blob_BBox PARAMS((BLOB *Blob));
  174.  
  175.  
  176.  
  177. /*****************************************************************************
  178. * Local variables
  179. ******************************************************************************/
  180.  
  181. METHODS Blob_Methods =
  182. {
  183.   All_Blob_Intersections,
  184.   Inside_Blob, Blob_Normal,
  185.   Copy_Blob,
  186.   Translate_Blob, Rotate_Blob, Scale_Blob, Transform_Blob,
  187.   Invert_Blob, Destroy_Blob
  188. };
  189.  
  190. static BSPHERE_TREE **Queue;
  191.  
  192. /* Maximum number of entries in queue during hierarchy traversal. */
  193.  
  194. static unsigned Max_Queue_Size = 1024;
  195.  
  196.  
  197.  
  198. /*****************************************************************************
  199. *
  200. * FUNCTION
  201. *
  202. *   All_Blob_Intersections
  203. *
  204. * INPUT
  205. *
  206. *   Object      - Object
  207. *   Ray         - Ray
  208. *
  209. * OUTPUT
  210. *
  211. *   Depth_Stack - Intersection stack
  212. *
  213. * RETURNS
  214. *
  215. *   int - TRUE, if a intersection was found
  216. *   
  217. * AUTHOR
  218. *
  219. *   Alexander Enzmann
  220. *   
  221. * DESCRIPTION
  222. *
  223. *   Generate intervals of influence for each component. After these
  224. *   are made, determine their aggregate effect on the ray. As the
  225. *   individual intervals are checked, a quartic is generated
  226. *   that represents the density at a particular point on the ray.
  227. *
  228. *   Explanation for spherical components:
  229. *
  230. *   After making the substitutions in MakeBlob, there is a formula
  231. *   for each component that has the form:
  232. *
  233. *      c0 * r^4 + c1 * r^2 + c2.
  234. *
  235. *   In order to determine the influence on the ray of all of the
  236. *   individual components, we start by determining the distance
  237. *   from any point on the ray to the specified point.  This can
  238. *   be found using the pythagorean theorem, using C as the center
  239. *   of this component, P as the start of the ray, and D as the
  240. *   direction of travel of the ray:
  241. *
  242. *      r^2 = (t * D + P - C) . (t * D + P - C)
  243. *
  244. *   we insert this equation for each appearance of r^2 in the
  245. *   components' formula, giving:
  246. *
  247. *      r^2 = D.D t^2 + 2 t D . (P - C) + (P - C) . (P - C)
  248. *
  249. *   Since the direction vector has been normalized, D.D = 1.
  250. *   Using the substitutions:
  251. *
  252. *      t0 = (P - C) . (P - C),
  253. *      t1 = D . (P - C)
  254. *
  255. *   We can write the formula as:
  256. *
  257. *      r^2 = t0 + 2 t t1 + t^2
  258. *
  259. *   Taking r^2 and substituting into the formula for this component
  260. *   of the Blob we get the formula:
  261. *
  262. *      density = c0 * (r^2)^2 + c1 * r^2 + c2,
  263. *
  264. *   or:
  265. *
  266. *      density = c0 * (t0 + 2 t t1 + t^2)^2 +
  267. *                c1 * (t0 + 2 t t1 + t^2) +
  268. *                c2
  269. *
  270. *   Expanding terms and collecting with respect to "t" gives:
  271. *
  272. *      t^4 * c0 +
  273. *      t^3 * 4 c0 t1 +
  274. *      t^2 * (c1 + 2 * c0 t0 + 4 c0 t1^2)
  275. *      t   * 2 (c1 t1 + 2 c0 t0 t1) +
  276. *            c2 + c1*t0 + c0*t0^2
  277. *
  278. *   This formula can now be solved for "t" by any of the quartic
  279. *   root solvers that are available.
  280. *
  281. * CHANGES
  282. *
  283. *   Jul 1994 : Added code for cylindrical and ellipsoidical blobs. [DB]
  284. *
  285. *   Oct 1994 : Added code to convert polynomial into a bezier curve for
  286. *              a quick test if an intersection exists in an interval. [DB]
  287. *
  288. *   Sep 1995 : Added code to avoid numerical problems with distant blobs. [DB]
  289. *
  290. ******************************************************************************/
  291.  
  292. static int All_Blob_Intersections(Object, Ray, Depth_Stack)
  293. OBJECT *Object;
  294. RAY *Ray;
  295. ISTACK *Depth_Stack;
  296. {
  297.   int i, j, cnt;
  298.   int root_count, in_flag;
  299.   int Intersection_Found = FALSE;
  300.   DBL t0, t1, t2, c0, c1, c2;
  301.   DBL dist, len, *fcoeffs, coeffs[5], roots[4];
  302.   DBL start_dist;
  303.   VECTOR P, D, V1, PP, DD;
  304.   VECTOR IPoint;
  305.   Blob_Interval *intervals;
  306.   BLOB_ELEMENT *Element;
  307.   BLOB *Blob = (BLOB *)Object;
  308.  
  309.   DBL l, m, w, newcoeffs[5], dk[5];
  310.  
  311.   Increase_Counter(stats[Ray_Blob_Tests]);
  312.  
  313.   /* Transform the ray into blob space. */
  314.  
  315.   if (Blob->Trans != NULL)
  316.   {
  317.     MInvTransPoint(P, Ray->Initial, Blob->Trans);
  318.     MInvTransDirection(D, Ray->Direction, Blob->Trans);
  319.  
  320.     VLength(len, D);
  321.     VInverseScaleEq(D, len);
  322.   }
  323.   else
  324.   {
  325.     Assign_Vector(P, Ray->Initial);
  326.     Assign_Vector(D, Ray->Direction);
  327.  
  328.     len = 1.0;
  329.   }
  330.  
  331.   /* Allocate temporary intersection list. */
  332.  
  333.   intervals = (Blob_Interval *)POV_MALLOC(2*Blob->Data->Number_Of_Components*sizeof(Blob_Interval), "temporary blob intervals");
  334.  
  335.   /* Get the intervals along the ray where each component has an effect. */
  336.   
  337.   if ((cnt = determine_influences(P, D, Blob, DEPTH_TOLERANCE, intervals)) == 0)
  338.   {
  339.     /* Ray doesn't hit any of the component elements. */
  340.     
  341.     POV_FREE(intervals);
  342.     
  343.     return (FALSE);
  344.   }
  345.  
  346.   /* To avoid numerical problems we start at the first interval. */
  347.  
  348.   if ((start_dist = intervals[0].bound) < Small_Tolerance)
  349.   {
  350.     start_dist = 0.0;
  351.   }
  352.  
  353.   for (i = 0; i < cnt; i++)
  354.   {
  355.     intervals[i].bound -= start_dist;
  356.   }
  357.  
  358.   /* Get the new starting point. */
  359.  
  360.   VAddScaledEq(P, start_dist, D);
  361.  
  362.   /* Clear out the coefficients. */
  363.  
  364.   coeffs[0] =
  365.   coeffs[1] =
  366.   coeffs[2] =
  367.   coeffs[3] = 0.0;
  368.   coeffs[4] = - Blob->Data->Threshold;
  369.   
  370.   /*
  371.    * Step through the list of intersection points, adding the
  372.    * influence of each component as it appears. 
  373.    */
  374.  
  375.   fcoeffs = NULL;
  376.   
  377.   for (i = in_flag = 0; i < cnt; i++)
  378.   {
  379.     if ((intervals[i].type & 1) == ENTERING)
  380.     {
  381.       /*
  382.        * Something is just starting to influence the ray, so calculate
  383.        * its coefficients and add them into the pot. 
  384.        */
  385.       
  386.       in_flag++;
  387.       
  388.       Element = intervals[i].Element;
  389.       
  390.       switch (Element->Type)
  391.       {
  392.         case BLOB_SPHERE:
  393.         
  394.           VSub(V1, P, Element->O);
  395.         
  396.           VDot(t0, V1, V1);
  397.           VDot(t1, V1, D);
  398.  
  399.           c0 = Element->c[0];
  400.           c1 = Element->c[1];
  401.           c2 = Element->c[2];
  402.         
  403.           fcoeffs = &(Element->f[0]);
  404.         
  405.           fcoeffs[0] = c0;
  406.           fcoeffs[1] = 4.0 * c0 * t1;
  407.           fcoeffs[2] = 2.0 * c0 * (2.0 * t1 * t1 + t0) + c1;
  408.           fcoeffs[3] = 2.0 * t1 * (2.0 * c0 * t0 + c1);
  409.           fcoeffs[4] = t0 * (c0 * t0 + c1) + c2;
  410.         
  411.           break;
  412.         
  413.         case BLOB_ELLIPSOID:
  414.         
  415.           MInvTransPoint(PP, P, Element->Trans);
  416.           MInvTransDirection(DD, D, Element->Trans);
  417.         
  418.           VSub(V1, PP, Element->O);
  419.         
  420.           VDot(t0, V1, V1);
  421.           VDot(t1, V1, DD);
  422.           VDot(t2, DD, DD);
  423.         
  424.           c0 = Element->c[0];
  425.           c1 = Element->c[1];
  426.           c2 = Element->c[2];
  427.         
  428.           fcoeffs = &(Element->f[0]);
  429.         
  430.           fcoeffs[0] = c0 * t2 * t2;
  431.           fcoeffs[1] = 4.0 * c0 * t1 * t2;
  432.           fcoeffs[2] = 2.0 * c0 * (2.0 * t1 * t1 + t0 * t2) + c1 * t2;
  433.           fcoeffs[3] = 2.0 * t1 * (2.0 * c0 * t0 + c1);
  434.           fcoeffs[4] = t0 * (c0 * t0 + c1) + c2;
  435.         
  436.           break;
  437.         
  438.         case BLOB_BASE_HEMISPHERE:
  439.         case BLOB_APEX_HEMISPHERE:
  440.         
  441.           MInvTransPoint(PP, P, Element->Trans);
  442.           MInvTransDirection(DD, D, Element->Trans);
  443.         
  444.           if (Element->Type == BLOB_APEX_HEMISPHERE)
  445.           {
  446.             PP[Z] -= Element->len;
  447.           }
  448.         
  449.           VDot(t0, PP, PP);
  450.           VDot(t1, PP, DD);
  451.           VDot(t2, DD, DD);
  452.  
  453.           c0 = Element->c[0];
  454.           c1 = Element->c[1];
  455.           c2 = Element->c[2];
  456.  
  457.           fcoeffs = &(Element->f[0]);
  458.  
  459.           fcoeffs[0] = c0 * t2 * t2;
  460.           fcoeffs[1] = 4.0 * c0 * t1 * t2;
  461.           fcoeffs[2] = 2.0 * c0 * (2.0 * t1 * t1 + t0 * t2) + c1 * t2;
  462.           fcoeffs[3] = 2.0 * t1 * (2.0 * c0 * t0 + c1);
  463.           fcoeffs[4] = t0 * (c0 * t0 + c1) + c2;
  464.         
  465.           break;
  466.         
  467.         case BLOB_CYLINDER:
  468.  
  469.           /* Transform ray into cylinder space. */
  470.         
  471.           MInvTransPoint(PP, P, Element->Trans);
  472.           MInvTransDirection(DD, D, Element->Trans);
  473.         
  474.           t0 = PP[X] * PP[X] + PP[Y] * PP[Y];
  475.           t1 = PP[X] * DD[X] + PP[Y] * DD[Y];
  476.           t2 = DD[X] * DD[X] + DD[Y] * DD[Y];
  477.         
  478.           c0 = Element->c[0];
  479.           c1 = Element->c[1];
  480.           c2 = Element->c[2];
  481.         
  482.           fcoeffs = &(Element->f[0]);
  483.           
  484.           fcoeffs[0] = c0 * t2 * t2;
  485.           fcoeffs[1] = 4.0 * c0 * t1 * t2;
  486.           fcoeffs[2] = 2.0 * c0 * (2.0 * t1 * t1 + t0 * t2) + c1 * t2;
  487.           fcoeffs[3] = 2.0 * t1 * (2.0 * c0 * t0 + c1);
  488.           fcoeffs[4] = t0 * (c0 * t0 + c1) + c2;
  489.  
  490.           break;
  491.         
  492.         default:
  493.         
  494.           Error("Unknown blob component in All_Blob_Intersections().\n");
  495.       }
  496.       
  497.       for (j = 0; j < 5; j++)
  498.       {
  499.         coeffs[j] += fcoeffs[j];
  500.       }
  501.     }
  502.     else
  503.     {
  504.       /* 
  505.        * We are losing the influence of a component -->
  506.        * subtract off its coefficients. 
  507.        */
  508.       
  509.       fcoeffs = intervals[i].Element->f;
  510.       
  511.       for (j = 0; j < 5; j++)
  512.       {
  513.         coeffs[j] -= fcoeffs[j];
  514.       }
  515.       
  516.       /* If no components are currently affecting the ray ---> skip ahead. */
  517.       
  518.       if (--in_flag == 0)
  519.       {
  520.         continue;
  521.       }
  522.     }
  523.     
  524.     /* 
  525.      * If the following intersection lies close to the current intersection
  526.      * then first add/subtract next region before testing. [DB 7/94] 
  527.      */
  528.     
  529.     if ((i + 1 < cnt) && (fabs(intervals[i].bound - intervals[i + 1].bound) < EPSILON))
  530.     {
  531.       continue;
  532.     }
  533.  
  534.     /*
  535.      * Transform polynomial in a way that the interval boundaries are moved
  536.      * to 0 and 1, i. e. the roots of interest are between 0 and 1. [DB 10/94]
  537.      */
  538.     
  539.     l = intervals[i].bound;
  540.     m = intervals[i+1].bound;
  541.  
  542.     w = m - l;
  543.  
  544.     newcoeffs[0] = coeffs[0] * w * w * w * w;
  545.     newcoeffs[1] = (coeffs[1] + 4.0 * coeffs[0] * l) * w * w * w;
  546.     newcoeffs[2] = (3.0 * l * (2.0 * coeffs[0] * l + coeffs[1]) + coeffs[2]) * w * w;
  547.     newcoeffs[3] = (2.0 * l * (l * (2.0 * coeffs[0] * l + 1.5 * coeffs[1]) + coeffs[2]) + coeffs[3]) * w;
  548.     newcoeffs[4] = l * (l * (l * (coeffs[0] * l + coeffs[1]) + coeffs[2]) + coeffs[3]) + coeffs[4];
  549.  
  550.     /* Calculate coefficients of corresponding bezier curve. [DB 10/94] */
  551.  
  552. /*
  553.     dk[0] = newcoeffs[4];
  554.     dk[1] = newcoeffs[4] + 0.25 * newcoeffs[3];
  555.     dk[2] = newcoeffs[4] + 0.50 * newcoeffs[3] + newcoeffs[2] / 6.0;
  556.     dk[3] = newcoeffs[4] + 0.75 * newcoeffs[3] + 0.5 * newcoeffs[2] + 0.25 * newcoeffs[1];
  557.     dk[4] = newcoeffs[4] + newcoeffs[3] + newcoeffs[2] + newcoeffs[1] + newcoeffs[0];
  558. */
  559.     dk[0] = newcoeffs[4];
  560.     dk[1] = newcoeffs[4] + (newcoeffs[3] / 4);
  561.     dk[2] = newcoeffs[4] + (newcoeffs[3] / 2) + newcoeffs[2] / 6.0;
  562.     dk[3] = newcoeffs[4] + 0.75 * newcoeffs[3] + (newcoeffs[2] / 2) + (newcoeffs[1] / 4);
  563.     dk[4] = newcoeffs[4] + newcoeffs[3] + newcoeffs[2] + newcoeffs[1] + newcoeffs[0];
  564.  
  565.     /*
  566.      * Skip this interval if the ray doesn't intersect the convex hull of the
  567.      * bezier curve, because no valid intersection will be found. [DB 10/94]
  568.      */
  569.  
  570.     if (((dk[0] >= 0.0) && (dk[1] >= 0.0) && (dk[2] >= 0.0) && (dk[3] >= 0.0) && (dk[4] >= 0.0)) ||
  571.         ((dk[0] <= 0.0) && (dk[1] <= 0.0) && (dk[2] <= 0.0) && (dk[3] <= 0.0) && (dk[4] <= 0.0)))
  572.     {
  573.       continue;
  574.     }
  575.  
  576.     /*
  577.      * Now we could do bezier clipping to find the roots
  578.      * but I have no idea how this works. [DB 2/95]
  579.      */
  580.  
  581.  
  582.     /* Solve polynomial. */
  583.  
  584.     root_count = Solve_Polynomial(4, coeffs, roots, Test_Flag(Blob, STURM_FLAG), ROOT_TOLERANCE);
  585.  
  586.     /* See if any of the roots are valid. */
  587.  
  588.     for (j = 0; j < root_count; j++)
  589.     {
  590.       dist = roots[j];
  591.  
  592.       /*
  593.        * First see if the root is in the interval of influence of
  594.        * the currently active components.
  595.        */
  596.  
  597.       if ((dist >= intervals[i].bound) && (dist <= intervals[i+1].bound) &&
  598.           (dist > DEPTH_TOLERANCE) && (dist < Max_Distance))
  599.       {
  600.         /* Correct distance. */
  601.  
  602.         dist += start_dist;
  603.  
  604.         dist /= len;
  605.  
  606.         VEvaluateRay(IPoint, Ray->Initial, dist, Ray->Direction);
  607.  
  608.         if (Point_In_Clip(IPoint, Object->Clip))
  609.         {
  610.           push_entry(dist, IPoint, Object, Depth_Stack);
  611.  
  612.           Intersection_Found = TRUE;
  613.         }
  614.       }
  615.     }
  616.  
  617.     /*
  618.      * If the blob isn't used inside a CSG and we have found at least
  619.      * one intersection then we are ready, because all possible intersections
  620.      * will be further away (we have a sorted list!). [DB 7/94]
  621.      */
  622.  
  623.     if (!(Blob->Type & IS_CHILD_OBJECT) && (Intersection_Found))
  624.     {
  625.       break;
  626.     }
  627.   }
  628.  
  629.   if (Intersection_Found)
  630.   {
  631.     Increase_Counter(stats[Ray_Blob_Tests_Succeeded]);
  632.   }
  633.  
  634.   /* Get rid off temporary intersection list. */
  635.  
  636.   POV_FREE(intervals);
  637.  
  638.   return (Intersection_Found);
  639. }
  640.  
  641.  
  642.  
  643. /*****************************************************************************
  644. *
  645. * FUNCTION
  646. *
  647. *   insert_hit
  648. *
  649. * INPUT
  650. *
  651. *   Blob      - Pointer to blob structure
  652. *   Element   - Element to insert
  653. *   t0, t1    - Intersection depths
  654. *
  655. * OUTPUT
  656. *
  657. *   intervals - Pointer to sorted list of hits
  658. *   cnt       - Number of hits in intervals
  659. *
  660. * RETURNS
  661. *
  662. * AUTHOR
  663. *
  664. *   Alexander Enzmann
  665. *
  666. * DESCRIPTION
  667. *
  668. *   -
  669. *
  670. * CHANGES
  671. *
  672. *   Oct 1994 : Modified to use memmove instead of loops for copying. [DB]
  673. *
  674. *   Sep 1995 : Changed to allow use of memcpy if memmove isn't available. [AED]
  675. *
  676. ******************************************************************************/
  677.  
  678. static void insert_hit(Element, t0, t1, intervals, cnt)
  679. BLOB_ELEMENT *Element;
  680. DBL t0, t1;
  681. Blob_Interval *intervals;
  682. int *cnt;
  683. {
  684.   int k;
  685.  
  686.   /*
  687.    * Store the points of intersection. Keep track of: whether this is
  688.    * the start or end point of the hit, which component was pierced
  689.    * by the ray, and the point along the ray that the hit occured at.
  690.    */
  691.  
  692.   for (k = 0; (k < *cnt) && (t0 > intervals[k].bound); k++);
  693.  
  694.   if (k < *cnt)
  695.   {
  696.     /*
  697.      * This hit point is smaller than one that already exists -->
  698.      * bump the rest and insert it here.
  699.      */
  700.  
  701. #ifdef NO_MEMMOVE
  702.     int memcnt;
  703.  
  704.     for (memcnt = *cnt; memcnt > k; memcnt--)
  705.     {
  706.       memcpy(&intervals[memcnt], &intervals[memcnt-1], sizeof(Blob_Interval));
  707.     }
  708. #else /* have MEMMOVE */
  709.     memmove(&intervals[k+1], &intervals[k], (*cnt-k)*sizeof(Blob_Interval));
  710. #endif
  711.  
  712.     /* We are entering the component. */
  713.  
  714.     intervals[k].type    = Element->Type | ENTERING;
  715.     intervals[k].bound   = t0;
  716.     intervals[k].Element = Element;
  717.  
  718.     (*cnt)++;
  719.  
  720.     for (k = k + 1; (k < *cnt) && (t1 > intervals[k].bound); k++);
  721.  
  722.     if (k < *cnt)
  723.     {
  724. #ifdef NO_MEMMOVE
  725.       for (memcnt = *cnt; memcnt > k; memcnt--)
  726.       {
  727.         memcpy(&intervals[memcnt], &intervals[memcnt-1], sizeof(Blob_Interval));
  728.       }
  729. #else /* have MEMMOVE */
  730.       memmove(&intervals[k+1], &intervals[k], (*cnt-k)*sizeof(Blob_Interval));
  731. #endif
  732.  
  733.       /* We are exiting the component. */
  734.  
  735.       intervals[k].type    = Element->Type | EXITING;
  736.       intervals[k].bound   = t1;
  737.       intervals[k].Element = Element;
  738.     }
  739.     else
  740.     {
  741.       /* We are exiting the component. */
  742.  
  743.       intervals[*cnt].type    = Element->Type | EXITING;
  744.       intervals[*cnt].bound   = t1;
  745.       intervals[*cnt].Element = Element;
  746.     }
  747.  
  748.     (*cnt)++;
  749.   }
  750.   else
  751.   {
  752.     /* Just plop the start and end points at the end of the list */
  753.  
  754.     /* We are entering the component. */
  755.  
  756.     intervals[*cnt].type    = Element->Type | ENTERING;
  757.     intervals[*cnt].bound   = t0;
  758.     intervals[*cnt].Element = Element;
  759.  
  760.     (*cnt)++;
  761.  
  762.     /* We are exiting the component. */
  763.  
  764.     intervals[*cnt].type    = Element->Type | EXITING;
  765.     intervals[*cnt].bound   = t1;
  766.     intervals[*cnt].Element = Element;
  767.  
  768.     (*cnt)++;
  769.   }
  770. }
  771.  
  772.  
  773.  
  774. /*****************************************************************************
  775. *
  776. * FUNCTION
  777. *
  778. *   intersect_cylinder
  779. *
  780. * INPUT
  781. *
  782. *   Element    - Pointer to element structure
  783. *   P, D       - Ray = P + t * D
  784. *   mindist    - Min. valid distance
  785. *
  786. * OUTPUT
  787. *
  788. *   tmin, tmax - Intersection depths found
  789. *
  790. * RETURNS
  791. *
  792. * AUTHOR
  793. *
  794. *   Dieter Bayer (with help from Alexander Enzmann)
  795. *
  796. * DESCRIPTION
  797. *
  798. *   -
  799. *
  800. * CHANGES
  801. *
  802. *   Jul 1994 : Creation.
  803. *
  804. ******************************************************************************/
  805.  
  806. static int intersect_cylinder(Element, P, D, mindist, tmin, tmax)
  807. BLOB_ELEMENT *Element;
  808. VECTOR P, D;
  809. DBL mindist, *tmin, *tmax;
  810. {
  811.   DBL a, b, c, d, t, u, v, w, len;
  812.   DBL rdiv;
  813.   VECTOR PP, DD;
  814.  
  815.   /* Transform ray into cylinder space. */
  816.  
  817.   MInvTransPoint(PP, P, Element->Trans);
  818.   MInvTransDirection(DD, D, Element->Trans);
  819.  
  820.   VLength(len, DD);
  821.   VInverseScaleEq(DD, len);
  822.  
  823.   /* Intersect ray with cylinder. */
  824.  
  825.   a = DD[X] * DD[X] + DD[Y] * DD[Y];
  826.  
  827.   if (a > EPSILON)
  828.   {
  829.     b = PP[X] * DD[X] + PP[Y] * DD[Y];
  830.     c = PP[X] * PP[X] + PP[Y] * PP[Y] - Element->rad2;
  831.  
  832.     d = b * b - a * c;
  833.  
  834.     if (d > EPSILON)
  835.     {
  836.       d = sqrt(d);
  837.  
  838.       rdiv = (1.0 / a);
  839.       t = ( - b + d) * rdiv;
  840. /*
  841.       t = ( - b + d) / a;
  842. */
  843.  
  844.       w = PP[Z] + t * DD[Z];
  845.  
  846.       if ((w >= 0.0) && (w <= Element->len))
  847.       {
  848.         if (t < *tmin) { *tmin = t; }
  849.         if (t > *tmax) { *tmax = t; }
  850.       }
  851.  
  852. /*
  853.       t = ( - b - d) / a;
  854. */
  855.       t = ( - b - d) * rdiv;
  856.  
  857.       w = PP[Z] + t * DD[Z];
  858.  
  859.       if ((w >= 0.0) && (w <= Element->len))
  860.       {
  861.         if (t < *tmin) { *tmin = t; }
  862.         if (t > *tmax) { *tmax = t; }
  863.       }
  864.     }
  865.   }
  866.  
  867.   /* Intersect base/cap plane. */
  868.  
  869.   if (fabs(DD[Z]) > EPSILON)
  870.   {
  871.     /* Intersect base plane. */
  872.  
  873.     rdiv = 1.0 / DD[Z];
  874.     t = - PP[Z] * rdiv;
  875. /*
  876.     t = - PP[Z] / DD[Z];
  877. */
  878.  
  879.     u = PP[X] + t * DD[X];
  880.     v = PP[Y] + t * DD[Y];
  881.  
  882.     if ((u * u + v * v) <= Element->rad2)
  883.     {
  884.       if (t < *tmin) { *tmin = t; }
  885.       if (t > *tmax) { *tmax = t; }
  886.     }
  887.  
  888.     /* Intersect cap plane. */
  889.  
  890. /*
  891.     t = (Element->len - PP[Z]) / DD[Z];
  892. */
  893.     t = (Element->len - PP[Z]) * rdiv;
  894.  
  895.     u = PP[X] + t * DD[X];
  896.     v = PP[Y] + t * DD[Y];
  897.  
  898.     if ((u * u + v * v) <= Element->rad2)
  899.     {
  900.       if (t < *tmin) { *tmin = t; }
  901.       if (t > *tmax) { *tmax = t; }
  902.     }
  903.   }
  904.  
  905.   /* Check if the intersections are valid. */
  906.  
  907.   rdiv = (1.0 / len);
  908.   *tmin *= rdiv;
  909.   *tmax *= rdiv;
  910.  
  911. /*
  912.   *tmin /= len;
  913.   *tmax /= len;
  914. */
  915.  
  916.   if (*tmin < mindist) { *tmin = 0.0; }
  917.   if (*tmax < mindist) { *tmax = 0.0; }
  918.  
  919.   if (*tmin >= *tmax)
  920.   {
  921.     return (FALSE);
  922.   }
  923.  
  924.   return (TRUE);
  925. }
  926.  
  927.  
  928.  
  929. /*****************************************************************************
  930. *
  931. * FUNCTION
  932. *
  933. *   intersect_ellipsoid
  934. *
  935. * INPUT
  936. *
  937. *   Element    - Pointer to element structure
  938. *   P, D       - Ray = P + t * D
  939. *   mindist    - Min. valid distance
  940. *
  941. * OUTPUT
  942. *
  943. *   tmin, tmax - Intersection depths found
  944. *
  945. * RETURNS
  946. *
  947. * AUTHOR
  948. *
  949. *   Dieter Bayer
  950. *
  951. * DESCRIPTION
  952. *
  953. *   -
  954. *
  955. * CHANGES
  956. *
  957. *   Sep 1994 : Creation.
  958. *
  959. ******************************************************************************/
  960.  
  961. static int intersect_ellipsoid(Element, P, D, mindist, tmin, tmax)
  962. BLOB_ELEMENT *Element;
  963. VECTOR P, D;
  964. DBL mindist, *tmin, *tmax;
  965. {
  966.   DBL b, d, t, len;
  967.   VECTOR V1, PP, DD;
  968.   DBL rdiv;
  969.  
  970.   MInvTransPoint(PP, P, Element->Trans);
  971.   MInvTransDirection(DD, D, Element->Trans);
  972.  
  973.   VLength(len, DD);
  974.   VInverseScaleEq(DD, len);
  975.  
  976.   VSub(V1, PP, Element->O);
  977.   VDot(b, V1, DD);
  978.   VDot(t, V1, V1);
  979.  
  980.   d = b * b - t + Element->rad2;
  981.  
  982.   if (d < EPSILON)
  983.   {
  984.     return (FALSE);
  985.   }
  986.  
  987.   d = sqrt(d);
  988.  
  989. /*
  990.   *tmax = ( - b + d) / len;  if (*tmax < mindist) { *tmax = 0.0; }
  991.   *tmin = ( - b - d) / len;  if (*tmin < mindist) { *tmin = 0.0; }
  992. */
  993.   rdiv = (1.0 / len);
  994.   *tmax = ( - b + d) * rdiv;  if (*tmax < mindist) { *tmax = 0.0; }
  995.   *tmin = ( - b - d) * rdiv;  if (*tmin < mindist) { *tmin = 0.0; }
  996.  
  997.   if (*tmax == *tmin)
  998.   {
  999.     return (FALSE);
  1000.   }
  1001.   else
  1002.   {
  1003.     if (*tmax < *tmin)
  1004.     {
  1005.       d = *tmin;  *tmin = *tmax;  *tmax = d;
  1006.     }
  1007.   }
  1008.  
  1009.   return (TRUE);
  1010. }
  1011.  
  1012.  
  1013.  
  1014. /*****************************************************************************
  1015. *
  1016. * FUNCTION
  1017. *
  1018. *   intersect_hemisphere
  1019. *
  1020. * INPUT
  1021. *
  1022. *   Element    - Pointer to element structure
  1023. *   P, D       - Ray = P + t * D
  1024. *   mindist    - Min. valid distance
  1025. *
  1026. * OUTPUT
  1027. *
  1028. *   tmin, tmax - Intersection depths found
  1029. *
  1030. * RETURNS
  1031. *
  1032. * AUTHOR
  1033. *
  1034. *   Dieter Bayer
  1035. *
  1036. * DESCRIPTION
  1037. *
  1038. *   -
  1039. *
  1040. * CHANGES
  1041. *
  1042. *   Jul 1994 : Creation (with help from Alexander Enzmann).
  1043. *
  1044. ******************************************************************************/
  1045.  
  1046. static int intersect_hemisphere(Element, P, D, mindist, tmin, tmax)
  1047. BLOB_ELEMENT *Element;
  1048. VECTOR P, D;
  1049. DBL mindist, *tmin, *tmax;
  1050. {
  1051.   DBL b, d, t, z1, z2, len;
  1052.   VECTOR PP, DD;
  1053.   DBL rdiv;
  1054.   
  1055.   /* Transform ray into hemisphere space. */
  1056.  
  1057.   MInvTransPoint(PP, P, Element->Trans);
  1058.   MInvTransDirection(DD, D, Element->Trans);
  1059.  
  1060.   VLength(len, DD);
  1061.   VInverseScaleEq(DD, len);
  1062.  
  1063.   if (Element->Type == BLOB_BASE_HEMISPHERE)
  1064.   {
  1065.     VDot(b, PP, DD);
  1066.     VDot(t, PP, PP);
  1067.  
  1068.     d = b * b - t + Element->rad2;
  1069.  
  1070.     if (d < EPSILON)
  1071.     {
  1072.       return (FALSE);
  1073.     }
  1074.  
  1075.     d = sqrt(d);
  1076.  
  1077.     *tmax = - b + d;
  1078.     *tmin = - b - d;
  1079.  
  1080.     if (*tmax < *tmin)
  1081.     {
  1082.       d = *tmin;  *tmin = *tmax;  *tmax = d;
  1083.     }
  1084.  
  1085.     /* Cut intersection at the plane. */
  1086.  
  1087.     z1 = PP[Z] + *tmin * DD[Z];
  1088.     z2 = PP[Z] + *tmax * DD[Z];
  1089.  
  1090.     /* If both points are inside --> no intersection */
  1091.  
  1092.     if ((z1 >= 0.0) && (z2 >= 0.0))
  1093.     {
  1094.       return (FALSE);
  1095.     }
  1096.  
  1097.     /* If both points are outside --> intersections found */
  1098.  
  1099.     if ((z1 < 0.0) && (z2 < 0.0))
  1100.     {
  1101.  
  1102. /*
  1103.       *tmin /= len;
  1104.       *tmax /= len;
  1105. */
  1106.       rdiv = (1.0 / len);
  1107.       *tmin *= rdiv;
  1108.       *tmax *= rdiv;
  1109.  
  1110.       return (TRUE);
  1111.     }
  1112.  
  1113.     /* Determine intersection with plane. */
  1114.  
  1115.     t = - PP[Z] / DD[Z];
  1116.  
  1117.     if (z1 >= 0.0)
  1118.     {
  1119.       /* Ray is crossing the plane from inside to outside. */
  1120.  
  1121.       *tmin = (t < mindist) ? 0.0 : t;
  1122.     }
  1123.     else
  1124.     {
  1125.       /* Ray is crossing the plane from outside to inside. */
  1126.  
  1127.       *tmax = (t < mindist) ? 0.0 : t;
  1128.     }
  1129.  
  1130.     rdiv = (1.0 / len);
  1131.     *tmin *= rdiv;
  1132.     *tmax *= rdiv;
  1133. /*
  1134.     *tmin /= len;
  1135.     *tmax /= len;
  1136. */
  1137.  
  1138.     return (TRUE);
  1139.   }
  1140.   else
  1141.   {
  1142.     PP[Z] -= Element->len;
  1143.  
  1144.     VDot(b, PP, DD);
  1145.     VDot(t, PP, PP);
  1146.  
  1147.     d = b * b - t + Element->rad2;
  1148.  
  1149.     if (d < EPSILON)
  1150.     {
  1151.       return (FALSE);
  1152.     }
  1153.  
  1154.     d = sqrt(d);
  1155.  
  1156.     *tmax = - b + d;
  1157.     *tmin = - b - d;
  1158.  
  1159.     if (*tmax < *tmin)
  1160.     {
  1161.       d = *tmin;  *tmin = *tmax;  *tmax = d;
  1162.     }
  1163.  
  1164.     /* Cut intersection at the plane. */
  1165.  
  1166.     z1 = PP[Z] + *tmin * DD[Z];
  1167.     z2 = PP[Z] + *tmax * DD[Z];
  1168.  
  1169.     /* If both points are inside --> no intersection */
  1170.  
  1171.     if ((z1 <= 0.0) && (z2 <= 0.0))
  1172.     {
  1173.       return (FALSE);
  1174.     }
  1175.  
  1176.     /* If both points are outside --> intersections found */
  1177.  
  1178.     if ((z1 > 0.0) && (z2 > 0.0))
  1179.     {
  1180.  
  1181.       rdiv = (1.0 / len);
  1182.       *tmin *= rdiv;
  1183.       *tmax *= rdiv;
  1184. /*
  1185.       *tmin /= len;
  1186.       *tmax /= len;
  1187. */
  1188.       return (TRUE);
  1189.     }
  1190.  
  1191.     /* Determine intersection with plane. */
  1192.  
  1193.     t = - PP[Z] / DD[Z];
  1194.  
  1195.     if (z1 <= 0.0)
  1196.     {
  1197.       /* Ray is crossing the plane from inside to outside. */
  1198.  
  1199.       *tmin = (t < mindist) ? 0.0 : t;
  1200.     }
  1201.     else
  1202.     {
  1203.       /* Ray is crossing the plane from outside to inside. */
  1204.  
  1205.       *tmax = (t < mindist) ? 0.0 : t;
  1206.     }
  1207.  
  1208.     rdiv = (1.0 / len);
  1209.     *tmin *= rdiv;
  1210.     *tmax *= rdiv;
  1211.  
  1212. /*
  1213.     *tmin /= len;
  1214.     *tmax /= len;
  1215. */
  1216.     return (TRUE);
  1217.   }
  1218. }
  1219.  
  1220.  
  1221.  
  1222. /*****************************************************************************
  1223. *
  1224. * FUNCTION
  1225. *
  1226. *   intersect_sphere
  1227. *
  1228. * INPUT
  1229. *
  1230. *   Element    - Pointer to element structure
  1231. *   P, D       - Ray = P + t * D
  1232. *   mindist    - Min. valid distance
  1233. *
  1234. * OUTPUT
  1235. *
  1236. *   tmin, tmax - Intersection depths found
  1237. *
  1238. * RETURNS
  1239. *
  1240. * AUTHOR
  1241. *
  1242. *   Dieter Bayer
  1243. *
  1244. * DESCRIPTION
  1245. *
  1246. *   -
  1247. *
  1248. * CHANGES
  1249. *
  1250. *   Jul 1994 : Creation (with help from Alexander Enzmann).
  1251. *
  1252. ******************************************************************************/
  1253.  
  1254. static int intersect_sphere(Element, P, D, mindist, tmin, tmax)
  1255. BLOB_ELEMENT *Element;
  1256. VECTOR P, D;
  1257. DBL mindist, *tmin, *tmax;
  1258. {
  1259.   DBL b, d, t;
  1260.   VECTOR V1;
  1261.  
  1262.   VSub(V1, P, Element->O);
  1263.   VDot(b, V1, D);
  1264.   VDot(t, V1, V1);
  1265.  
  1266.   d = b * b - t + Element->rad2;
  1267.  
  1268.   if (d < EPSILON)
  1269.   {
  1270.     return (FALSE);
  1271.   }
  1272.  
  1273.   d = sqrt(d);
  1274.  
  1275.   *tmax = - b + d;  if (*tmax < mindist) { *tmax = 0.0; }
  1276.   *tmin = - b - d;  if (*tmin < mindist) { *tmin = 0.0; }
  1277.  
  1278.   if (*tmax == *tmin)
  1279.   {
  1280.     return (FALSE);
  1281.   }
  1282.   else
  1283.   {
  1284.     if (*tmax < *tmin)
  1285.     {
  1286.       d = *tmin;  *tmin = *tmax;  *tmax = d;
  1287.     }
  1288.   }
  1289.  
  1290.   return (TRUE);
  1291. }
  1292.  
  1293.  
  1294.  
  1295. /*****************************************************************************
  1296. *
  1297. * FUNCTION
  1298. *
  1299. *   intersect_element
  1300. *
  1301. * INPUT
  1302. *
  1303. *   P, D       - Ray = P + t * D
  1304. *   Element    - Pointer to element structure
  1305. *   mindist    - Min. valid distance
  1306. *
  1307. * OUTPUT
  1308. *
  1309. *   tmin, tmax - Intersection depths found
  1310. *
  1311. * RETURNS
  1312. *
  1313. * AUTHOR
  1314. *
  1315. *   Dieter Bayer
  1316. *
  1317. * DESCRIPTION
  1318. *
  1319. *   -
  1320. *
  1321. * CHANGES
  1322. *
  1323. *   Jul 1994 : Creation.
  1324. *
  1325. ******************************************************************************/
  1326.  
  1327. static int intersect_element(P, D, Element, mindist, tmin, tmax)
  1328. VECTOR P, D;
  1329. BLOB_ELEMENT *Element;
  1330. DBL mindist, *tmin, *tmax;
  1331. {
  1332. #ifdef BLOB_EXTRA_STATS
  1333.   Increase_Counter(stats[Blob_Element_Tests]);
  1334. #endif
  1335.  
  1336.   *tmin = BOUND_HUGE;
  1337.   *tmax = - BOUND_HUGE;
  1338.  
  1339.   switch (Element->Type)
  1340.   {
  1341.     case BLOB_SPHERE:
  1342.  
  1343.       if (!intersect_sphere(Element, P, D, mindist, tmin, tmax))
  1344.       {
  1345.        return (FALSE);
  1346.       }
  1347.  
  1348.       break;
  1349.  
  1350.     case BLOB_ELLIPSOID:
  1351.  
  1352.       if (!intersect_ellipsoid(Element, P, D, mindist, tmin, tmax))
  1353.       {
  1354.         return (FALSE);
  1355.       }
  1356.  
  1357.       break;
  1358.  
  1359.     case BLOB_BASE_HEMISPHERE:
  1360.     case BLOB_APEX_HEMISPHERE:
  1361.  
  1362.       if (!intersect_hemisphere(Element, P, D, mindist, tmin, tmax))
  1363.       {
  1364.         return (FALSE);
  1365.       }
  1366.  
  1367.       break;
  1368.  
  1369.     case BLOB_CYLINDER:
  1370.  
  1371.       if (!intersect_cylinder(Element, P, D, mindist, tmin, tmax))
  1372.       {
  1373.         return (FALSE);
  1374.       }
  1375.  
  1376.       break;
  1377.   }
  1378.  
  1379. #ifdef BLOB_EXTRA_STATS
  1380.   Increase_Counter(stats[Blob_Element_Tests_Succeeded]);
  1381. #endif
  1382.  
  1383.   return (TRUE);
  1384. }
  1385.  
  1386.  
  1387.  
  1388. /*****************************************************************************
  1389. *
  1390. * FUNCTION
  1391. *
  1392. *   determine_influences
  1393. *
  1394. * INPUT
  1395. *
  1396. *   P, D       - Ray = P + t * D
  1397. *   Blob       - Pointer to blob structure
  1398. *   mindist    - Min. valid distance
  1399. *
  1400. * OUTPUT
  1401. *
  1402. *   intervals  - Sorted list of intersections found
  1403. *
  1404. * RETURNS
  1405. *
  1406. *   int - Number of intersection found
  1407. *
  1408. * AUTHOR
  1409. *
  1410. *   Dieter Bayer
  1411. *
  1412. * DESCRIPTION
  1413. *
  1414. *   Make a sorted list of points along the ray at which the various blob
  1415. *   components start and stop adding their influence.
  1416. *
  1417. * CHANGES
  1418. *
  1419. *   Jul 1994 : Added code for bounding hierarchy traversal. [DB]
  1420. *
  1421. ******************************************************************************/
  1422.  
  1423. static int determine_influences(P, D, Blob, mindist, intervals)
  1424. VECTOR P, D;
  1425. BLOB *Blob;
  1426. DBL mindist;
  1427. Blob_Interval *intervals;
  1428. {
  1429.   int i;
  1430.   int cnt, size;
  1431.   DBL b, t, t0, t1;
  1432.   VECTOR V1;
  1433.   BSPHERE_TREE *Tree;
  1434.  
  1435.   cnt = 0;
  1436.  
  1437.   if (Blob->Data->Tree == NULL)
  1438.   {
  1439.     /* There's no bounding hierarchy so just step through all elements. */
  1440.  
  1441.     for (i = 0; i < Blob->Data->Number_Of_Components; i++)
  1442.     {
  1443.       if (intersect_element(P, D, &Blob->Data->Entry[i], mindist, &t0, &t1))
  1444.       {
  1445.         insert_hit(&Blob->Data->Entry[i], t0, t1, intervals, &cnt);
  1446.       }
  1447.     }
  1448.   }
  1449.   else
  1450.   {
  1451.     /* Use blob's bounding hierarchy. */
  1452.  
  1453.     size = 0;
  1454.  
  1455.     Queue[size++] = Blob->Data->Tree;
  1456.  
  1457.     while (size > 0)
  1458.     {
  1459.       Tree = Queue[--size];
  1460.  
  1461.       /* Test if current node is a leaf. */
  1462.  
  1463.       if (Tree->Entries <= 0)
  1464.       {
  1465.         /* Test element. */
  1466.  
  1467.         if (intersect_element(P, D, (BLOB_ELEMENT *)Tree->Node, mindist, &t0, &t1))
  1468.         {
  1469.           insert_hit((BLOB_ELEMENT *)Tree->Node, t0, t1, intervals, &cnt);
  1470.         }
  1471.       }
  1472.       else
  1473.       {
  1474.         /* Test all sub-nodes. */
  1475.  
  1476.         for (i = 0; i < (int)Tree->Entries; i++)
  1477.         {
  1478. #ifdef BLOB_EXTRA_STATS
  1479.           Increase_Counter(stats[Blob_Bound_Tests]);
  1480. #endif
  1481.  
  1482.           VSub(V1, Tree->Node[i]->C, P);
  1483.           VDot(b, V1, D);
  1484.           VDot(t, V1, V1);
  1485.  
  1486.           if ((t - Sqr(b)) <= Tree->Node[i]->r2)
  1487.           {
  1488. #ifdef BLOB_EXTRA_STATS
  1489.             Increase_Counter(stats[Blob_Bound_Tests_Succeeded]);
  1490. #endif
  1491.  
  1492.             insert_node(Tree->Node[i], &size);
  1493.           }
  1494.         }
  1495.       }
  1496.     }
  1497.   }
  1498.  
  1499.   return (cnt);
  1500. }
  1501.  
  1502.  
  1503.  
  1504. /*****************************************************************************
  1505. *
  1506. * FUNCTION
  1507. *
  1508. *   calculate_element_field
  1509. *
  1510. * INPUT
  1511. *
  1512. *   Element - Pointer to element structure
  1513. *   P       - Point whos field value is calculated
  1514. *
  1515. * OUTPUT
  1516. *
  1517. * RETURNS
  1518. *
  1519. *   DBL - Field value
  1520. *
  1521. * AUTHOR
  1522. *
  1523. *   Alexander Enzmann
  1524. *
  1525. * DESCRIPTION
  1526. *
  1527. *   Calculate the field value of a single element in a given point P
  1528. *   (which must already have been transformed into blob space).
  1529. *
  1530. * CHANGES
  1531. *
  1532. *   Jul 1994 : Added code for cylindrical and ellipsoidical blobs. [DB]
  1533. *
  1534. ******************************************************************************/
  1535.  
  1536. static DBL calculate_element_field(Element, P)
  1537. BLOB_ELEMENT *Element;
  1538. VECTOR P;
  1539. {
  1540.   DBL rad2, density;
  1541.   VECTOR V1, PP;
  1542.  
  1543.   density = 0.0;
  1544.  
  1545.   switch (Element->Type)
  1546.   {
  1547.     case BLOB_SPHERE:
  1548.  
  1549.       VSub(V1, P, Element->O);
  1550.  
  1551.       VDot(rad2, V1, V1);
  1552.  
  1553.       if (rad2 < Element->rad2)
  1554.       {
  1555.         density = rad2 * (rad2 * Element->c[0] + Element->c[1]) + Element->c[2];
  1556.       }
  1557.  
  1558.       break;
  1559.  
  1560.     case BLOB_ELLIPSOID:
  1561.  
  1562.       MInvTransPoint(PP, P, Element->Trans);
  1563.  
  1564.       VSub(V1, PP, Element->O);
  1565.  
  1566.       VDot(rad2, V1, V1);
  1567.  
  1568.       if (rad2 < Element->rad2)
  1569.       {
  1570.         density = rad2 * (rad2 * Element->c[0] + Element->c[1]) + Element->c[2];
  1571.       }
  1572.  
  1573.       break;
  1574.  
  1575.     case BLOB_BASE_HEMISPHERE:
  1576.  
  1577.       MInvTransPoint(PP, P, Element->Trans);
  1578.  
  1579.       if (PP[Z] <= 0.0)
  1580.       {
  1581.         VDot(rad2, PP, PP);
  1582.  
  1583.         if (rad2 <= Element->rad2)
  1584.         {
  1585.           density = rad2 * (rad2 * Element->c[0] + Element->c[1]) + Element->c[2];
  1586.         }
  1587.       }
  1588.  
  1589.       break;
  1590.  
  1591.     case BLOB_APEX_HEMISPHERE:
  1592.  
  1593.       MInvTransPoint(PP, P, Element->Trans);
  1594.  
  1595.       PP[Z] -= Element->len;
  1596.  
  1597.       if (PP[Z] >= 0.0)
  1598.       {
  1599.         VDot(rad2, PP, PP);
  1600.  
  1601.         if (rad2 <= Element->rad2)
  1602.         {
  1603.           density = rad2 * (rad2 * Element->c[0] + Element->c[1]) + Element->c[2];
  1604.         }
  1605.       }
  1606.  
  1607.       break;
  1608.  
  1609.     case BLOB_CYLINDER:
  1610.  
  1611.       MInvTransPoint(PP, P, Element->Trans);
  1612.  
  1613.       if ((PP[Z] >= 0.0) && (PP[Z] <= Element->len))
  1614.       {
  1615.         if ((rad2 = Sqr(PP[X]) + Sqr(PP[Y])) <= Element->rad2)
  1616.         {
  1617.           density = rad2 * (rad2 * Element->c[0] + Element->c[1]) + Element->c[2];
  1618.         }
  1619.       }
  1620.  
  1621.       break;
  1622.   }
  1623.  
  1624.   return (density);
  1625. }
  1626.  
  1627.  
  1628.  
  1629. /*****************************************************************************
  1630. *
  1631. * FUNCTION
  1632. *
  1633. *   calculate_field_value
  1634. *
  1635. * INPUT
  1636. *
  1637. *   Blob - Pointer to blob structure
  1638. *   P       - Point whos field value is calculated
  1639. *
  1640. * OUTPUT
  1641. *
  1642. * RETURNS
  1643. *
  1644. *   DBL - Field value
  1645. *
  1646. * AUTHOR
  1647. *
  1648. *   Dieter Bayer
  1649. *
  1650. * DESCRIPTION
  1651. *
  1652. *   Calculate the field value of a blob in a given point P
  1653. *   (which must already have been transformed into blob space).
  1654. *
  1655. * CHANGES
  1656. *
  1657. *   Jul 1994 : Added code for bounding hierarchy traversal. [DB]
  1658. *
  1659. ******************************************************************************/
  1660.  
  1661. static DBL calculate_field_value(Blob, P)
  1662. BLOB *Blob;
  1663. VECTOR P;
  1664. {
  1665.   int i;
  1666.   int size;
  1667.   DBL density, rad2;
  1668.   VECTOR V1;
  1669.   BSPHERE_TREE *Tree;
  1670.  
  1671.   density = 0.0;
  1672.  
  1673.   if (Blob->Data->Tree == NULL)
  1674.   {
  1675.     /* There's no tree --> step through all elements. */
  1676.  
  1677.     for (i = 0; i < Blob->Data->Number_Of_Components; i++)
  1678.     {
  1679.       density += calculate_element_field(&Blob->Data->Entry[i], P);
  1680.     }
  1681.   }
  1682.   else
  1683.   {
  1684.     /* A tree exists --> step through the tree. */
  1685.  
  1686.     size = 0;
  1687.  
  1688.     Queue[size++] = Blob->Data->Tree;
  1689.  
  1690.     while (size > 0)
  1691.     {
  1692.       Tree = Queue[--size];
  1693.  
  1694.       /* Test if current node is a leaf. */
  1695.  
  1696.       if (Tree->Entries <= 0)
  1697.       {
  1698.         density += calculate_element_field((BLOB_ELEMENT *)Tree->Node, P);
  1699.       }
  1700.       else
  1701.       {
  1702.         /* Test all sub-nodes. */
  1703.  
  1704.         for (i = 0; i < (int)Tree->Entries; i++)
  1705.         {
  1706.           /* Insert sub-node if we are inside. */
  1707.  
  1708.           VSub(V1, P, Tree->Node[i]->C);
  1709.  
  1710.           VDot(rad2, V1, V1);
  1711.  
  1712.           if (rad2 <= Tree->Node[i]->r2)
  1713.           {
  1714.             insert_node(Tree->Node[i], &size);
  1715.           }
  1716.         }
  1717.       }
  1718.     }
  1719.   }
  1720.  
  1721.   return (density);
  1722. }
  1723.  
  1724.  
  1725.  
  1726. /*****************************************************************************
  1727. *
  1728. * FUNCTION
  1729. *
  1730. *   Inside_Blob
  1731. *
  1732. * INPUT
  1733. *
  1734. *   Test_Point - Point to test
  1735. *   Object     - Pointer to blob structure
  1736. *
  1737. * OUTPUT
  1738. *
  1739. * RETURNS
  1740. *
  1741. *   int - TRUE if Test_Point is inside
  1742. *
  1743. * AUTHOR
  1744. *
  1745. *   Alexander Enzmann
  1746. *
  1747. * DESCRIPTION
  1748. *
  1749. *   Calculate the density at the given point and then compare to
  1750. *   the threshold to see if we are in or out of the blob.
  1751. *
  1752. * CHANGES
  1753. *
  1754. *   -
  1755. *
  1756. ******************************************************************************/
  1757.  
  1758. static int Inside_Blob(Test_Point, Object)
  1759. VECTOR Test_Point;
  1760. OBJECT *Object;
  1761. {
  1762.   VECTOR New_Point;
  1763.   BLOB *Blob = (BLOB *) Object;
  1764.  
  1765.   /* Transform the point into blob space. */
  1766.  
  1767.   if (Blob->Trans != NULL)
  1768.   {
  1769.     MInvTransPoint(New_Point, Test_Point, Blob->Trans);
  1770.   }
  1771.   else
  1772.   {
  1773.     Assign_Vector(New_Point, Test_Point);
  1774.   }
  1775.  
  1776.   if (calculate_field_value(Blob, New_Point) > Blob->Data->Threshold - INSIDE_TOLERANCE)
  1777.   {
  1778.     /* We are inside. */
  1779.  
  1780.     return (!Test_Flag(Blob, INVERTED_FLAG));
  1781.   }
  1782.   else
  1783.   {
  1784.     /* We are outside. */
  1785.  
  1786.     return (Test_Flag(Blob, INVERTED_FLAG));
  1787.   }
  1788. }
  1789.  
  1790.  
  1791.  
  1792. /*****************************************************************************
  1793. *
  1794. * FUNCTION
  1795. *
  1796. *   element_normal
  1797. *
  1798. * INPUT
  1799. *
  1800. *   P       - Surface point
  1801. *   Element - Pointer to element structure
  1802. *
  1803. * OUTPUT
  1804. *
  1805. *   Result  - Element's normal
  1806. *
  1807. * RETURNS
  1808. *
  1809. * AUTHOR
  1810. *
  1811. *   Dieter Bayer
  1812. *
  1813. * DESCRIPTION
  1814. *
  1815. *   Calculate the normal of a single element in the point P.
  1816. *
  1817. * CHANGES
  1818. *
  1819. *   Jul 1994 : Creation (with help from Alexander Enzmann).
  1820. *
  1821. ******************************************************************************/
  1822.  
  1823. static void element_normal(Result, P, Element)
  1824. VECTOR Result, P;
  1825. BLOB_ELEMENT *Element;
  1826. {
  1827.   DBL val, dist;
  1828.   VECTOR V1, PP;
  1829.  
  1830.   switch (Element->Type)
  1831.   {
  1832.     case BLOB_SPHERE:
  1833.  
  1834.       VSub(V1, P, Element->O);
  1835.  
  1836.       VDot(dist, V1, V1);
  1837.  
  1838.       if (dist <= Element->rad2)
  1839.       {
  1840.         val = -2.0 * Element->c[0] * dist - Element->c[1];
  1841.  
  1842.         VAddScaledEq(Result, val, V1);
  1843.       }
  1844.  
  1845.       break;
  1846.  
  1847.     case BLOB_ELLIPSOID:
  1848.  
  1849.       MInvTransPoint(PP, P, Element->Trans);
  1850.  
  1851.       VSub(V1, PP, Element->O);
  1852.  
  1853.       VDot(dist, V1, V1);
  1854.  
  1855.       if (dist <= Element->rad2)
  1856.       {
  1857.         val = -2.0 * Element->c[0] * dist - Element->c[1];
  1858.  
  1859.         MTransNormal(V1, V1, Element->Trans);
  1860.  
  1861.         VAddScaledEq(Result, val, V1);
  1862.       }
  1863.  
  1864.       break;
  1865.  
  1866.     case BLOB_BASE_HEMISPHERE:
  1867.  
  1868.       MInvTransPoint(PP, P, Element->Trans);
  1869.  
  1870.       if (PP[Z] <= 0.0)
  1871.       {
  1872.         VDot(dist, PP, PP);
  1873.  
  1874.         if (dist <= Element->rad2)
  1875.         {
  1876.           val = -2.0 * Element->c[0] * dist - Element->c[1];
  1877.  
  1878.           MTransNormal(PP, PP, Element->Trans);
  1879.  
  1880.           VAddScaledEq(Result, val, PP);
  1881.         }
  1882.       }
  1883.  
  1884.       break;
  1885.  
  1886.     case BLOB_APEX_HEMISPHERE:
  1887.  
  1888.       MInvTransPoint(PP, P, Element->Trans);
  1889.  
  1890.       PP[Z] -= Element->len;
  1891.  
  1892.       if (PP[Z] >= 0.0)
  1893.       {
  1894.         VDot(dist, PP, PP);
  1895.  
  1896.         if (dist <= Element->rad2)
  1897.         {
  1898.           val = -2.0 * Element->c[0] * dist - Element->c[1];
  1899.  
  1900.           MTransNormal(PP, PP, Element->Trans);
  1901.  
  1902.           VAddScaledEq(Result, val, PP);
  1903.         }
  1904.       }
  1905.  
  1906.       break;
  1907.  
  1908.     case BLOB_CYLINDER:
  1909.  
  1910.       MInvTransPoint(PP, P, Element->Trans);
  1911.  
  1912.       if ((PP[Z] >= 0.0) && (PP[Z] <= Element->len))
  1913.       {
  1914.         if ((dist = Sqr(PP[X]) + Sqr(PP[Y])) <= Element->rad2)
  1915.         {
  1916.           val = -2.0 * Element->c[0] * dist - Element->c[1];
  1917.  
  1918.           PP[Z] = 0.0;
  1919.  
  1920.           MTransNormal(PP, PP, Element->Trans);
  1921.  
  1922.           VAddScaledEq(Result, val, PP);
  1923.         }
  1924.       }
  1925.  
  1926.       break;
  1927.   }
  1928. }
  1929.  
  1930.  
  1931.  
  1932. /*****************************************************************************
  1933. *
  1934. * FUNCTION
  1935. *
  1936. *   Blob_Normal
  1937. *
  1938. * INPUT
  1939. *
  1940. *   Object  - Pointer to blob structure
  1941. *   Inter   - Pointer to intersection
  1942. *
  1943. * OUTPUT
  1944. *
  1945. *   Result  - Blob's normal
  1946. *
  1947. * RETURNS
  1948. *
  1949. * AUTHOR
  1950. *
  1951. *   Alexander Enzmann
  1952. *
  1953. * DESCRIPTION
  1954. *
  1955. *   Calculate the blob's surface normal in the intersection point.
  1956. *
  1957. * CHANGES
  1958. *
  1959. *   Jul 1994 : Added code for bounding hierarchy traversal. [DB]
  1960. *
  1961. ******************************************************************************/
  1962.  
  1963. static void Blob_Normal(Result, Object, Inter)
  1964. OBJECT *Object;
  1965. VECTOR Result;
  1966. INTERSECTION *Inter;
  1967. {
  1968.   int i;
  1969.   int size;
  1970.   DBL dist, val;
  1971.   VECTOR New_Point, V1;
  1972.   BLOB *Blob = (BLOB *) Object;
  1973.   BSPHERE_TREE *Tree;
  1974.  
  1975.   /* Transform the point into the blob space. */
  1976.  
  1977.   if (Blob->Trans != NULL)
  1978.   {
  1979.     MInvTransPoint(New_Point, Inter->IPoint, Blob->Trans);
  1980.   }
  1981.   else
  1982.   {
  1983.     Assign_Vector(New_Point, Inter->IPoint);
  1984.   }
  1985.  
  1986.   Make_Vector(Result, 0.0, 0.0, 0.0);
  1987.  
  1988.   /* For each component that contributes to this point, add its bit to the normal */
  1989.  
  1990.   if (Blob->Data->Tree == NULL)
  1991.   {
  1992.     /* There's no tree --> step through all elements. */
  1993.  
  1994.     for (i = 0; i < Blob->Data->Number_Of_Components; i++)
  1995.     {
  1996.       element_normal(Result, New_Point, &(Blob->Data->Entry[i]));
  1997.     }
  1998.   }
  1999.   else
  2000.   {
  2001.     /* A tree exists --> step through the tree. */
  2002.  
  2003.     size = 0;
  2004.  
  2005.     Queue[size++] = Blob->Data->Tree;
  2006.  
  2007.     while (size > 0)
  2008.     {
  2009.       Tree = Queue[--size];
  2010.  
  2011.       /* Test if current node is a leaf. */
  2012.  
  2013.       if (Tree->Entries <= 0)
  2014.       {
  2015.         element_normal(Result, New_Point, (BLOB_ELEMENT *)Tree->Node);
  2016.       }
  2017.       else
  2018.       {
  2019.         /* Test all sub-nodes. */
  2020.  
  2021.         for (i = 0; i < (int)Tree->Entries; i++)
  2022.         {
  2023.           /* Insert sub-node if we are inside. */
  2024.  
  2025.           VSub(V1, New_Point, Tree->Node[i]->C);
  2026.  
  2027.           VDot(dist, V1, V1);
  2028.  
  2029.           if (dist <= Tree->Node[i]->r2)
  2030.           {
  2031.             insert_node(Tree->Node[i], &size);
  2032.           }
  2033.         }
  2034.       }
  2035.     }
  2036.   }
  2037.  
  2038.   VDot(val, Result, Result);
  2039.  
  2040.   if (val == 0.0)
  2041.   {
  2042.     Make_Vector(Result, 1.0, 0.0, 0.0);
  2043.   }
  2044.   else
  2045.   {
  2046.     /* Normalize normal vector. */
  2047.  
  2048.     val = 1.0 / sqrt(val);
  2049.  
  2050.     VScaleEq(Result, val);
  2051.   }
  2052.  
  2053.   /* Transform back to world space. */
  2054.  
  2055.   if (Blob->Trans != NULL)
  2056.   {
  2057.     MTransNormal(Result, Result, Blob->Trans);
  2058.  
  2059.     VNormalize(Result, Result);
  2060.   }
  2061. }
  2062.  
  2063.  
  2064.  
  2065. /*****************************************************************************
  2066. *
  2067. * FUNCTION
  2068. *
  2069. *   Translate_Blob
  2070. *
  2071. * INPUT
  2072. *
  2073. *   Vector - Translation vector
  2074. *
  2075. * OUTPUT
  2076. *
  2077. *   Object - Pointer to blob structure
  2078. *
  2079. * RETURNS
  2080. *
  2081. * AUTHOR
  2082. *
  2083. *   Alexander Enzmann
  2084. *
  2085. * DESCRIPTION
  2086. *
  2087. *   Translate a blob.
  2088. *
  2089. * CHANGES
  2090. *
  2091. *   -
  2092. *
  2093. ******************************************************************************/
  2094.  
  2095. static void Translate_Blob(Object, Vector, Trans)
  2096. OBJECT *Object;
  2097. VECTOR Vector;
  2098. TRANSFORM *Trans;
  2099. {
  2100.   Transform_Blob(Object, Trans);
  2101. }
  2102.  
  2103.  
  2104.  
  2105. /*****************************************************************************
  2106. *
  2107. * FUNCTION
  2108. *
  2109. *   Rotate_Blob
  2110. *
  2111. * INPUT
  2112. *
  2113. *   Vector - Rotation vector
  2114. *
  2115. * OUTPUT
  2116. *
  2117. *   Object - Pointer to blob structure
  2118. *
  2119. * RETURNS
  2120. *
  2121. * AUTHOR
  2122. *
  2123. *   Alexander Enzmann
  2124. *
  2125. * DESCRIPTION
  2126. *
  2127. *   Rotate a blob.
  2128. *
  2129. * CHANGES
  2130. *
  2131. *   -
  2132. *
  2133. ******************************************************************************/
  2134.  
  2135. static void Rotate_Blob(Object, Vector, Trans)
  2136. OBJECT *Object;
  2137. VECTOR Vector;
  2138. TRANSFORM *Trans;
  2139. {
  2140.   Transform_Blob(Object, Trans);
  2141. }
  2142.  
  2143.  
  2144.  
  2145. /*****************************************************************************
  2146. *
  2147. * FUNCTION
  2148. *
  2149. *   Scale_Blob
  2150. *
  2151. * INPUT
  2152. *
  2153. *   Vector - Scaling vector
  2154. *
  2155. * OUTPUT
  2156. *
  2157. *   Object - Pointer to blob structure
  2158. *
  2159. * RETURNS
  2160. *
  2161. * AUTHOR
  2162. *
  2163. *   Alexander Enzmann
  2164. *
  2165. * DESCRIPTION
  2166. *
  2167. *   Scale a blob.
  2168. *
  2169. * CHANGES
  2170. *
  2171. *   -
  2172. *
  2173. ******************************************************************************/
  2174.  
  2175. static void Scale_Blob(Object, Vector, Trans)
  2176. OBJECT *Object;
  2177. VECTOR Vector;
  2178. TRANSFORM *Trans;
  2179. {
  2180.   Transform_Blob(Object, Trans);
  2181. }
  2182.  
  2183.  
  2184.  
  2185. /*****************************************************************************
  2186. *
  2187. * FUNCTION
  2188. *
  2189. *   Transform_Blob
  2190. *
  2191. * INPUT
  2192. *
  2193. *   Trans  - Pointer to transformation
  2194. *
  2195. * OUTPUT
  2196. *
  2197. *   Object - Pointer to blob structure
  2198. *
  2199. * RETURNS
  2200. *   
  2201. * AUTHOR
  2202. *
  2203. *   Alexander Enzmann
  2204. *   
  2205. * DESCRIPTION
  2206. *
  2207. *   Transform a blob.
  2208. *
  2209. * CHANGES
  2210. *
  2211. *   -
  2212. *
  2213. ******************************************************************************/
  2214.  
  2215. static void Transform_Blob(Object, Trans)
  2216. OBJECT *Object;
  2217. TRANSFORM *Trans;
  2218. {
  2219.   int i;
  2220.   BLOB *Blob = (BLOB *)Object;
  2221.  
  2222.   if (Blob->Trans == NULL)
  2223.   {
  2224.     Blob->Trans = Create_Transform();
  2225.   }
  2226.  
  2227.   Recompute_BBox(&Object->BBox, Trans);
  2228.  
  2229.   Compose_Transforms(Blob->Trans, Trans);
  2230.  
  2231.   for (i = 0; i < Blob->Data->Number_Of_Components; i++)
  2232.   {
  2233.     Transform_Textures(Blob->Element_Texture[i], Trans);
  2234.   }
  2235. }
  2236.  
  2237.  
  2238.  
  2239. /*****************************************************************************
  2240. *
  2241. * FUNCTION
  2242. *
  2243. *   Invert_Blob
  2244. *
  2245. * INPUT
  2246. *
  2247. *   Object - Pointer to blob structure
  2248. *
  2249. * OUTPUT
  2250. *
  2251. *   Object
  2252. *
  2253. * RETURNS
  2254. *
  2255. * AUTHOR
  2256. *
  2257. *   Alexander Enzmann
  2258. *
  2259. * DESCRIPTION
  2260. *
  2261. *   Invert a blob.
  2262. *
  2263. * CHANGES
  2264. *
  2265. *   -
  2266. *
  2267. ******************************************************************************/
  2268.  
  2269. static void Invert_Blob(Object)
  2270. OBJECT *Object;
  2271. {
  2272.   Invert_Flag(Object, INVERTED_FLAG);
  2273. }
  2274.  
  2275.  
  2276.  
  2277. /*****************************************************************************
  2278. *
  2279. * FUNCTION
  2280. *
  2281. *   Create_Blob
  2282. *
  2283. * INPUT
  2284. *
  2285. *   Object - Pointer to blob structure
  2286. *
  2287. * OUTPUT
  2288. *
  2289. *   Object
  2290. *
  2291. * RETURNS
  2292. *
  2293. * AUTHOR
  2294. *
  2295. *   Alexander Enzmann
  2296. *
  2297. * DESCRIPTION
  2298. *
  2299. *   Create a new blob.
  2300. *
  2301. * CHANGES
  2302. *
  2303. *   -
  2304. *
  2305. ******************************************************************************/
  2306.  
  2307. BLOB *Create_Blob()
  2308. {
  2309.   BLOB *New;
  2310.  
  2311.   New = (BLOB *)POV_MALLOC(sizeof (BLOB), "blob");
  2312.  
  2313.   INIT_OBJECT_FIELDS(New, BLOB_OBJECT, &Blob_Methods)
  2314.  
  2315.   Set_Flag(New, HIERARCHY_FLAG);
  2316.  
  2317.   New->Trans = NULL;
  2318.  
  2319.   New->Data = NULL;
  2320.  
  2321.   New->Element_Texture = NULL;
  2322.  
  2323.   return (New);
  2324. }
  2325.  
  2326.  
  2327.  
  2328. /*****************************************************************************
  2329. *
  2330. * FUNCTION
  2331. *
  2332. *   Copy_Blob
  2333. *
  2334. * INPUT
  2335. *
  2336. *   Object - Pointer to blob structure
  2337. *
  2338. * OUTPUT
  2339. *
  2340. *   Object
  2341. *
  2342. * RETURNS
  2343. *
  2344. * AUTHOR
  2345. *
  2346. *   Alexander Enzmann
  2347. *
  2348. * DESCRIPTION
  2349. *
  2350. *   Copy a blob.
  2351. *
  2352. *   NOTE: The components are not copied, only the number of references is
  2353. *         counted, so that Destroy_Blob() knows if they can be destroyed.
  2354. *
  2355. * CHANGES
  2356. *
  2357. *   Jul 1994 : Added code for blob data reference counting. [DB]
  2358. *
  2359. ******************************************************************************/
  2360.  
  2361. static void *Copy_Blob(Object)
  2362. OBJECT *Object;
  2363. {
  2364.   int i;
  2365.   BLOB *New, *Old = (BLOB *)Object;
  2366.  
  2367.   New = Create_Blob();
  2368.  
  2369.   /* Copy blob. */
  2370.  
  2371.   *New = *Old;
  2372.  
  2373.   New->Trans = Copy_Transform(New->Trans);
  2374.  
  2375.   New->Data->References++;
  2376.  
  2377.   New->Element_Texture = (TEXTURE **)POV_MALLOC(New->Data->Number_Of_Components*sizeof(TEXTURE *), "blob texture list");
  2378.  
  2379.   for (i = 0; i < New->Data->Number_Of_Components; i++)
  2380.   {
  2381.     New->Element_Texture[i] = Copy_Textures(Old->Element_Texture[i]);
  2382.   }
  2383.  
  2384.   return (New);
  2385. }
  2386.  
  2387.  
  2388.  
  2389. /*****************************************************************************
  2390. *
  2391. * FUNCTION
  2392. *
  2393. *   Create_Blob_List_Element
  2394. *
  2395. * INPUT
  2396. *
  2397. * OUTPUT
  2398. *
  2399. * RETURNS
  2400. *
  2401. *   blobstackptr - Pointer to blob element
  2402. *
  2403. * AUTHOR
  2404. *
  2405. *   Dieter Bayer
  2406. *
  2407. * DESCRIPTION
  2408. *
  2409. *   Create a new blob element in the component list used during parsing.
  2410. *
  2411. * CHANGES
  2412. *
  2413. *   Sep 1994 : Creation.
  2414. *
  2415. ******************************************************************************/
  2416.  
  2417. blobstackptr Create_Blob_List_Element()
  2418. {
  2419.   blobstackptr New;
  2420.  
  2421.   New = (blobstackptr)POV_MALLOC(sizeof(struct blob_list_struct), "blob component");
  2422.  
  2423.   init_blob_element(&New->elem);
  2424.  
  2425.   return (New);
  2426. }
  2427.  
  2428.  
  2429.  
  2430. /*****************************************************************************
  2431. *
  2432. * FUNCTION
  2433. *
  2434. *   Destroy_Blob
  2435. *
  2436. * INPUT
  2437. *
  2438. *   Object - Pointer to blob structure
  2439. *
  2440. * OUTPUT
  2441. *
  2442. *   Object
  2443. *
  2444. * RETURNS
  2445. *
  2446. * AUTHOR
  2447. *
  2448. *   Alexander Enzmann
  2449. *
  2450. * DESCRIPTION
  2451. *
  2452. *   Destroy a blob.
  2453. *
  2454. *   NOTE: The blob data is destroyed if they are no longer used by any copy.
  2455. *
  2456. * CHANGES
  2457. *
  2458. *   Jul 1994 : Added code for blob data reference counting. [DB]
  2459. *
  2460. *   Dec 1994 : Fixed memory leakage. [DB]
  2461. *
  2462. *   Aug 1995 : Fixed freeing of already freed memory. [DB]
  2463. *
  2464. ******************************************************************************/
  2465.  
  2466. static void Destroy_Blob(Object)
  2467. OBJECT *Object;
  2468. {
  2469.   int i;
  2470.   BLOB *Blob = (BLOB *)Object;
  2471.  
  2472.   Destroy_Transform(Blob->Trans);
  2473.  
  2474.   for (i = 0; i < Blob->Data->Number_Of_Components; i++)
  2475.   {
  2476.     Destroy_Textures(Blob->Element_Texture[i]);
  2477.   }
  2478.  
  2479.   POV_FREE(Blob->Element_Texture);
  2480.  
  2481.   if (--(Blob->Data->References) == 0)
  2482.   {
  2483.     Destroy_Bounding_Sphere_Hierarchy(Blob->Data->Tree);
  2484.  
  2485.     for (i = 0; i < Blob->Data->Number_Of_Components; i++)
  2486.     {
  2487.       /*
  2488.        * Make sure to destroy multiple references of a texture
  2489.        * and/or transformation only once. Multiple references
  2490.        * are only used with cylindrical blobs. Thus it's
  2491.        * enough to ignore all cylinder caps.
  2492.        */
  2493.  
  2494.       if ((Blob->Data->Entry[i].Type == BLOB_SPHERE) ||
  2495.           (Blob->Data->Entry[i].Type == BLOB_ELLIPSOID) ||
  2496.           (Blob->Data->Entry[i].Type == BLOB_CYLINDER))
  2497.       {
  2498.         Destroy_Transform(Blob->Data->Entry[i].Trans);
  2499.  
  2500.         Blob->Data->Entry[i].Trans   = NULL;
  2501.         Blob->Data->Entry[i].Texture = NULL;
  2502.       }
  2503.     }
  2504.  
  2505.     POV_FREE(Blob->Data->Entry);
  2506.  
  2507.     POV_FREE(Blob->Data);
  2508.   }
  2509.  
  2510.   POV_FREE(Object);
  2511. }
  2512.  
  2513.  
  2514.  
  2515. /*****************************************************************************
  2516. *
  2517. * FUNCTION
  2518. *
  2519. *   Compute_Blob_BBox
  2520. *
  2521. * INPUT
  2522. *
  2523. *   Blob - Blob
  2524. *
  2525. * OUTPUT
  2526. *
  2527. *   Blob
  2528. *
  2529. * RETURNS
  2530. *
  2531. * AUTHOR
  2532. *
  2533. *   Dieter Bayer
  2534. *
  2535. * DESCRIPTION
  2536. *
  2537. *   Calculate the bounding box of a blob.
  2538. *
  2539. * CHANGES
  2540. *
  2541. *   Aug 1994 : Creation.
  2542. *
  2543. ******************************************************************************/
  2544.  
  2545. static void Compute_Blob_BBox(Blob)
  2546. BLOB *Blob;
  2547. {
  2548.   int i;
  2549.   DBL radius, radius2;
  2550.   VECTOR Center, Min, Max;
  2551.  
  2552.   Make_Vector(Min, BOUND_HUGE, BOUND_HUGE, BOUND_HUGE);
  2553.   Make_Vector(Max, - BOUND_HUGE, - BOUND_HUGE, - BOUND_HUGE);
  2554.  
  2555.   for (i = 0; i < Blob->Data->Number_Of_Components; i++)
  2556.   {
  2557.     get_element_bounding_sphere(&Blob->Data->Entry[i], Center, &radius2);
  2558.  
  2559.     radius = sqrt(radius2);
  2560.  
  2561.     Min[X] = min(Min[X], Center[X] - radius);
  2562.     Min[Y] = min(Min[Y], Center[Y] - radius);
  2563.     Min[Z] = min(Min[Z], Center[Z] - radius);
  2564.     Max[X] = max(Max[X], Center[X] + radius);
  2565.     Max[Y] = max(Max[Y], Center[Y] + radius);
  2566.     Max[Z] = max(Max[Z], Center[Z] + radius);
  2567.   }
  2568.  
  2569.   Make_BBox_from_min_max(Blob->BBox, Min, Max);
  2570.  
  2571.   if (Blob->Trans != NULL)
  2572.   {
  2573.     Recompute_BBox(&Blob->BBox, Blob->Trans);
  2574.   }
  2575. }
  2576.  
  2577.  
  2578.  
  2579. /*****************************************************************************
  2580. *
  2581. * FUNCTION
  2582. *
  2583. *   get_element_bounding_sphere
  2584. *
  2585. * INPUT
  2586. *
  2587. *   Element - Pointer to element
  2588. *   Center  - Bounding sphere's center
  2589. *   Radius2 - Bounding sphere's squared radius
  2590. *
  2591. * OUTPUT
  2592. *
  2593. *   Center, Radius2
  2594. *
  2595. * RETURNS
  2596. *
  2597. * AUTHOR
  2598. *
  2599. *   Dieter Bayer
  2600. *
  2601. * DESCRIPTION
  2602. *
  2603. *   Calculate the bounding sphere of a blob element.
  2604. *
  2605. * CHANGES
  2606. *
  2607. *   Sep 1994 : Creation.
  2608. *
  2609. ******************************************************************************/
  2610.  
  2611. static void get_element_bounding_sphere(Element, Center, Radius2)
  2612. BLOB_ELEMENT *Element;
  2613. VECTOR Center;
  2614. DBL *Radius2;
  2615. {
  2616.   DBL r, r2 = 0.0;
  2617.   VECTOR C, H;
  2618.  
  2619.   switch (Element->Type)
  2620.   {
  2621.     case BLOB_SPHERE:
  2622.     case BLOB_ELLIPSOID:
  2623.  
  2624.       r2 = Element->rad2;
  2625.  
  2626.       Assign_Vector(C, Element->O);
  2627.  
  2628.       break;
  2629.  
  2630.     case BLOB_BASE_HEMISPHERE:
  2631.  
  2632.       r2 = Element->rad2;
  2633.  
  2634.       Make_Vector(C, 0.0, 0.0, 0.0);
  2635.  
  2636.       break;
  2637.  
  2638.     case BLOB_APEX_HEMISPHERE:
  2639.  
  2640.       r2 = Element->rad2;
  2641.  
  2642.       Make_Vector(C, 0.0, 0.0, Element->len);
  2643.  
  2644.       break;
  2645.  
  2646.     case BLOB_CYLINDER :
  2647.  
  2648.       Make_Vector(C, 0.0, 0.0, 0.5 * Element->len);
  2649.  
  2650.       r2 = Element->rad2 + Sqr(0.5 * Element->len);
  2651.  
  2652.       break;
  2653.   }
  2654.  
  2655.   /* Transform bounding sphere if necessary. */
  2656.  
  2657.   if (Element->Trans != NULL)
  2658.   {
  2659.     r = sqrt(r2);
  2660.  
  2661.     MTransPoint(C, C, Element->Trans);
  2662.  
  2663.     Make_Vector(H, r, r, r);
  2664.  
  2665.     MTransDirection(H, H, Element->Trans);
  2666.  
  2667.     r = max(max(fabs(H[X]), fabs(H[Y])), fabs(H[Z]));
  2668.  
  2669.     r2 = Sqr(r) + EPSILON;
  2670.   }
  2671.  
  2672.   Assign_Vector(Center, C);
  2673.  
  2674.   *Radius2 = r2;
  2675. }
  2676.  
  2677.  
  2678.  
  2679. /*****************************************************************************
  2680. *
  2681. * FUNCTION
  2682. *
  2683. *   init_blob_element
  2684. *
  2685. * INPUT
  2686. *
  2687. *   Element - Pointer to blob element
  2688. *
  2689. * OUTPUT
  2690. *
  2691. *   Element
  2692. *
  2693. * RETURNS
  2694. *
  2695. * AUTHOR
  2696. *
  2697. *   Dieter Bayer
  2698. *
  2699. * DESCRIPTION
  2700. *
  2701. *   Init blob element.
  2702. *
  2703. * CHANGES
  2704. *
  2705. *   Sep 1994 : Creation.
  2706. *
  2707. ******************************************************************************/
  2708.  
  2709. static void init_blob_element(Element)
  2710. BLOB_ELEMENT *Element;
  2711. {
  2712.   Element->Type = 0;
  2713.  
  2714.   Element->index = 0;
  2715.  
  2716.   Element->len =
  2717.   Element->rad2 = 0.0;
  2718.  
  2719.   Element->c[0] =
  2720.   Element->c[1] =
  2721.   Element->c[2] =
  2722.   Element->f[0] =
  2723.   Element->f[1] =
  2724.   Element->f[2] =
  2725.   Element->f[3] =
  2726.   Element->f[4] = 0.0;
  2727.  
  2728.   Make_Vector(Element->O, 0.0, 0.0, 0.0);
  2729.  
  2730.   Element->Texture = NULL;
  2731.  
  2732.   Element->Trans = NULL;
  2733. }
  2734.  
  2735.  
  2736.  
  2737. /*****************************************************************************
  2738. *
  2739. * FUNCTION
  2740. *
  2741. *   Make_Blob
  2742. *
  2743. * INPUT
  2744. *
  2745. *   Blob       - Pointer to blob structure
  2746. *   threshold  - Blob's threshold
  2747. *   BlobList   - Pointer to elements
  2748. *   npoints    - Number of elements
  2749. *
  2750. * OUTPUT
  2751. *
  2752. *   Blob
  2753. *
  2754. * RETURNS
  2755. *
  2756. * AUTHOR
  2757. *
  2758. *   Alexander Enzmann
  2759. *
  2760. * DESCRIPTION
  2761. *
  2762. *   Create a blob after it was read from the scene file.
  2763. *
  2764. *   Starting with the density function: (1-r^2)^2, we have a field
  2765. *   that varies in strength from 1 at r = 0 to 0 at r = 1. By
  2766. *   substituting r/rad for r, we can adjust the range of influence
  2767. *   of a particular component. By multiplication by coeff, we can
  2768. *   adjust the amount of total contribution, giving the formula:
  2769. *
  2770. *     coeff * (1 - (r/rad)^2)^2
  2771. *
  2772. *   This varies in strength from coeff at r = 0, to 0 at r = rad.
  2773. *
  2774. * CHANGES
  2775. *
  2776. *   Jul 1994 : Added code for cylindrical and ellipsoidical blobs. [DB]
  2777. *
  2778. ******************************************************************************/
  2779.  
  2780. void Make_Blob(Blob, threshold, BlobList, npoints)
  2781. BLOB *Blob;
  2782. DBL threshold;
  2783. blobstackptr BlobList;
  2784. int npoints;
  2785. {
  2786.   int i, count;
  2787.   DBL rad2, coeff;
  2788.   blobstackptr temp;
  2789.   BLOB_ELEMENT *Entry;
  2790.  
  2791.   if (npoints < 1)
  2792.   {
  2793.     Error("Need at least one component in a blob.");
  2794.   }
  2795.  
  2796.   /* Figure out how many components there will be. */
  2797.  
  2798.   temp = BlobList;
  2799.  
  2800.   for (i = count = 0; i < npoints; i++)
  2801.   {
  2802.     if (temp->elem.Type & BLOB_CYLINDER)
  2803.     {
  2804.       count += 3;
  2805.     }
  2806.     else
  2807.     {
  2808.       count++;
  2809.     }
  2810.  
  2811.     temp = temp->next;
  2812.  
  2813.     /* Test for too many components. [DB 12/94] */
  2814.  
  2815.     if (count >= MAX_BLOB_COMPONENTS)
  2816.     {
  2817.       Error("There are more than %d components in a blob.\n", MAX_BLOB_COMPONENTS);
  2818.     }
  2819.   }
  2820.  
  2821.   /* Initialize the blob data. */
  2822.  
  2823.   Blob->Data->Threshold = threshold;
  2824.  
  2825.   Entry = Blob->Data->Entry;
  2826.  
  2827.   for (i = 0; i < npoints; i++)
  2828.   {
  2829.     temp = BlobList;
  2830.  
  2831.     if ((fabs(temp->elem.c[2]) < EPSILON) || (temp->elem.rad2 < EPSILON))
  2832.     {
  2833.       Warning(0.0, "Degenerate Blob element\n");
  2834.     }
  2835.  
  2836.     /* Initialize component. */
  2837.  
  2838.     *Entry = temp->elem;
  2839.  
  2840.     /* We have a multi-texture blob. */
  2841.  
  2842.     if (Entry->Texture != NULL)
  2843.     {
  2844.       Set_Flag(Blob, MULTITEXTURE_FLAG);
  2845.     }
  2846.  
  2847.     /* Store blob specific information. */
  2848.  
  2849.     switch (temp->elem.Type)
  2850.     {
  2851.       case BLOB_ELLIPSOID :
  2852.       case BLOB_SPHERE :
  2853.  
  2854.         rad2 = temp->elem.rad2;
  2855.  
  2856.         coeff = temp->elem.c[2];
  2857.  
  2858.         Entry->c[0] = coeff / (rad2 * rad2);
  2859.         Entry->c[1] = -(2.0 * coeff) / rad2;
  2860.         Entry->c[2] = coeff;
  2861.  
  2862.         Entry++;
  2863.  
  2864.         break;
  2865.  
  2866.       case BLOB_CYLINDER :
  2867.  
  2868.         rad2 = temp->elem.rad2;
  2869.  
  2870.         coeff = temp->elem.c[2];
  2871.  
  2872.         /* Create cylindrical component. */
  2873.  
  2874.         Entry->c[0] = coeff / (rad2 * rad2);
  2875.         Entry->c[1] = -(2.0 * coeff) / rad2;
  2876.         Entry->c[2] = coeff;
  2877.  
  2878.         Entry++;
  2879.  
  2880.         /* Create hemispherical component at the base. */
  2881.  
  2882.         *Entry = temp->elem;
  2883.  
  2884.         Entry->Type = BLOB_BASE_HEMISPHERE;
  2885.  
  2886.         Entry->c[0] = coeff / (rad2 * rad2);
  2887.         Entry->c[1] = -(2.0 * coeff) / rad2;
  2888.         Entry->c[2] = coeff;
  2889.  
  2890.         Entry++;
  2891.  
  2892.         /* Create hemispherical component at the apex. */
  2893.  
  2894.         *Entry = temp->elem;
  2895.  
  2896.         Entry->Type = BLOB_APEX_HEMISPHERE;
  2897.  
  2898.         Entry->c[0] = coeff / (rad2 * rad2);
  2899.         Entry->c[1] = -(2.0 * coeff) / rad2;
  2900.         Entry->c[2] = coeff;
  2901.  
  2902.         Entry++;
  2903.  
  2904.         break;
  2905.  
  2906.       default :
  2907.  
  2908.         Error("Unknown blob component.\n");
  2909.     }
  2910.  
  2911.     /* Get rid of texture non longer needed. */
  2912.  
  2913.     BlobList = BlobList->next;
  2914.  
  2915.     Destroy_Textures(temp->elem.Texture);
  2916.  
  2917.     POV_FREE(temp);
  2918.   }
  2919.  
  2920.   for (i = 0; i < count; i++)
  2921.   {
  2922.     Blob->Data->Entry[i].index = i;
  2923.   }
  2924.  
  2925.   /* Compute bounding box. */
  2926.  
  2927.   Compute_Blob_BBox(Blob);
  2928.  
  2929.   /* Create bounding sphere hierarchy. */
  2930.  
  2931.   if (Test_Flag(Blob, HIERARCHY_FLAG))
  2932.   {
  2933.     build_bounding_hierarchy(Blob);
  2934.   }
  2935. }
  2936.  
  2937.  
  2938.  
  2939. /*****************************************************************************
  2940. *
  2941. * FUNCTION
  2942. *
  2943. *   Test_Blob_Opacity
  2944. *
  2945. * INPUT
  2946. *
  2947. *   Blob - Pointer to blob structure
  2948. *
  2949. * OUTPUT
  2950. *
  2951. *   Blob
  2952. *
  2953. * RETURNS
  2954. *
  2955. * AUTHOR
  2956. *
  2957. *   Dieter Bayer
  2958. *
  2959. * DESCRIPTION
  2960. *
  2961. *   Set the opacity flag of the blob according to the opacity
  2962. *   of the blob's texture(s).
  2963. *
  2964. * CHANGES
  2965. *
  2966. *   Apr 1996 : Creation.
  2967. *
  2968. ******************************************************************************/
  2969.  
  2970. void Test_Blob_Opacity(Blob)
  2971. BLOB *Blob;
  2972. {
  2973.   int i;
  2974.  
  2975.   /* Initialize opacity flag to the opacity of the object's texture. */
  2976.  
  2977.   if ((Blob->Texture == NULL) || (Test_Opacity(Blob->Texture)))
  2978.   {
  2979.     Set_Flag(Blob, OPAQUE_FLAG);
  2980.   }
  2981.  
  2982.   if (Test_Flag(Blob, MULTITEXTURE_FLAG))
  2983.   {
  2984.     for (i = 0; i < Blob->Data->Number_Of_Components; i++)
  2985.     {
  2986.       if (Blob->Element_Texture[i] != NULL)
  2987.       {
  2988.         /* If component's texture isn't opaque the blob is neither. */
  2989.  
  2990.         if (!Test_Opacity(Blob->Element_Texture[i]))
  2991.         {
  2992.           Clear_Flag(Blob, OPAQUE_FLAG);
  2993.         }
  2994.       }
  2995.     }
  2996.   }
  2997. }
  2998.  
  2999.  
  3000.  
  3001. /*****************************************************************************
  3002. *
  3003. * FUNCTION
  3004. *
  3005. *   build_bounding_hierarchy
  3006. *
  3007. * INPUT
  3008. *
  3009. * OUTPUT
  3010. *
  3011. * RETURNS
  3012. *
  3013. * AUTHOR
  3014. *
  3015. *   Dieter Bayer
  3016. *
  3017. * DESCRIPTION
  3018. *
  3019. *   Create the bounding sphere hierarchy.
  3020. *
  3021. * CHANGES
  3022. *
  3023. *   Oct 1994 : Creation. (Derived from the bounding slab creation code)
  3024. *
  3025. ******************************************************************************/
  3026.  
  3027. static void build_bounding_hierarchy(Blob)
  3028. BLOB *Blob;
  3029. {
  3030.   int i, nElem, maxelements;
  3031.   BSPHERE_TREE **Elements;
  3032.  
  3033.   nElem = (int)Blob->Data->Number_Of_Components;
  3034.  
  3035.   maxelements = 2 * nElem;
  3036.  
  3037.   /*
  3038.    * Now allocate an array to hold references to these elements.
  3039.    */
  3040.  
  3041.   Elements = (BSPHERE_TREE **)POV_MALLOC(maxelements*sizeof(BSPHERE_TREE *), "blob bounding hierarchy");
  3042.  
  3043.   /* Init list with blob elements. */
  3044.  
  3045.   for (i = 0; i < nElem; i++)
  3046.   {
  3047.     Elements[i] = (BSPHERE_TREE *)POV_MALLOC(sizeof(BSPHERE_TREE), "blob bounding hierarchy");
  3048.  
  3049.     Elements[i]->Entries = 0;
  3050.     Elements[i]->Node    = (BSPHERE_TREE **)&Blob->Data->Entry[i];
  3051.  
  3052.     get_element_bounding_sphere(&Blob->Data->Entry[i], Elements[i]->C, &Elements[i]->r2);
  3053.   }
  3054.  
  3055.   Build_Bounding_Sphere_Hierarchy(&Blob->Data->Tree, nElem, Elements);
  3056.  
  3057.   /* Get rid of the Elements array. */
  3058.  
  3059.   POV_FREE(Elements);
  3060. }
  3061.  
  3062.  
  3063.  
  3064. /*****************************************************************************
  3065. *
  3066. * FUNCTION
  3067. *
  3068. *   Determine_Blob_Textures
  3069. *
  3070. * INPUT
  3071. *
  3072. * OUTPUT
  3073. *
  3074. * RETURNS
  3075. *
  3076. * AUTHOR
  3077. *
  3078. *   Dieter Bayer
  3079. *
  3080. * DESCRIPTION
  3081. *
  3082. *   Determine the textures and weights of all components affecting
  3083. *   the given intersection point. The weights are calculated from
  3084. *   the field values and sum to 1.
  3085. *
  3086. * CHANGES
  3087. *
  3088. *   Sep 1994 : Creation.
  3089. *
  3090. *   Mar 1996 : Make the call to resize the textures/weights list just once
  3091. *              at the beginning instead of doing it for every element. [DB]
  3092. *
  3093. ******************************************************************************/
  3094.  
  3095. void Determine_Blob_Textures(Blob, IPoint, Count, Textures, Weights)
  3096. BLOB *Blob;
  3097. VECTOR IPoint;
  3098. int *Count;
  3099. TEXTURE **Textures;
  3100. DBL *Weights;
  3101. {
  3102.   int i;
  3103.   int size;
  3104.   DBL rad2, sum;
  3105.   VECTOR V1, P;
  3106.   BLOB_ELEMENT *Element;
  3107.   BSPHERE_TREE *Tree;
  3108.  
  3109.   /* Make sure we have enough room in the textures/weights list. */
  3110.  
  3111.   Reinitialize_Lighting_Code(Blob->Data->Number_Of_Components, &Textures, &Weights);
  3112.  
  3113.   /* Transform the point into the blob space. */
  3114.  
  3115.   if (Blob->Trans != NULL)
  3116.   {
  3117.     MInvTransPoint(P, IPoint, Blob->Trans);
  3118.   }
  3119.   else
  3120.   {
  3121.     Assign_Vector(P, IPoint);
  3122.   }
  3123.  
  3124.   *Count = 0;
  3125.  
  3126.   if (Blob->Data->Tree == NULL)
  3127.   {
  3128.     /* There's no tree --> step through all elements. */
  3129.  
  3130.     for (i = 0; i < Blob->Data->Number_Of_Components; i++)
  3131.     {
  3132.       Element = &Blob->Data->Entry[i];
  3133.  
  3134.       determine_element_texture(Blob, Element, Blob->Element_Texture[i], P, Count, Textures, Weights);
  3135.     }
  3136.   }
  3137.   else
  3138.   {
  3139.     /* A tree exists --> step through the tree. */
  3140.  
  3141.     size = 0;
  3142.  
  3143.     Queue[size++] = Blob->Data->Tree;
  3144.  
  3145.     while (size > 0)
  3146.     {
  3147.       Tree = Queue[--size];
  3148.  
  3149.       /* Test if current node is a leaf. */
  3150.  
  3151.       if (Tree->Entries <= 0)
  3152.       {
  3153.         determine_element_texture(Blob, (BLOB_ELEMENT *)Tree->Node, Blob->Element_Texture[((BLOB_ELEMENT *)Tree->Node)->index], P, Count, Textures, Weights);
  3154.       }
  3155.       else
  3156.       {
  3157.         /* Test all sub-nodes. */
  3158.  
  3159.         for (i = 0; i < (int)Tree->Entries; i++)
  3160.         {
  3161.           /* Insert sub-node if we are inside. */
  3162.  
  3163.           VSub(V1, P, Tree->Node[i]->C);
  3164.  
  3165.           VDot(rad2, V1, V1);
  3166.  
  3167.           if (rad2 <= Tree->Node[i]->r2)
  3168.           {
  3169.             insert_node(Tree->Node[i], &size);
  3170.           }
  3171.         }
  3172.       }
  3173.     }
  3174.   }
  3175.  
  3176.   /* Normalize weights so that their sum is 1. */
  3177.  
  3178.   if (*Count > 0)
  3179.   {
  3180.     sum = 0.0;
  3181.  
  3182.     for (i = 0; i < *Count; i++)
  3183.     {
  3184.       sum += Weights[i];
  3185.     }
  3186.  
  3187.     if (sum > 0.0)
  3188.     {
  3189.       for (i = 0; i < *Count; i++)
  3190.       {
  3191.         Weights[i] /= sum;
  3192.       }
  3193.     }
  3194.   }
  3195. }
  3196.  
  3197.  
  3198.  
  3199. /*****************************************************************************
  3200. *
  3201. * FUNCTION
  3202. *
  3203. *   determine_element_texture
  3204. *
  3205. * INPUT
  3206. *
  3207. * OUTPUT
  3208. *
  3209. * RETURNS
  3210. *
  3211. * AUTHOR
  3212. *
  3213. *   Dieter Bayer
  3214. *
  3215. * DESCRIPTION
  3216. *
  3217. *   If the intersection point is inside the component calculate
  3218. *   the field density and store the element's texture and the field
  3219. *   value in the texture/weight list.
  3220. *
  3221. * CHANGES
  3222. *
  3223. *   Sep 1994 : Creation.
  3224. *
  3225. ******************************************************************************/
  3226.  
  3227. static void determine_element_texture(Blob, Element, Texture, P, Count, Textures, Weights)
  3228. BLOB *Blob;
  3229. BLOB_ELEMENT *Element;
  3230. VECTOR P;
  3231. int *Count;
  3232. TEXTURE *Texture, **Textures;
  3233. DBL *Weights;
  3234. {
  3235.   int i;
  3236.   DBL density;
  3237.  
  3238.   density = fabs(calculate_element_field(Element, P));
  3239.  
  3240.   if (density > 0.0)
  3241.   {
  3242.     if (Texture == NULL)
  3243.     {
  3244.       Textures[*Count] = Blob->Texture;
  3245.     }
  3246.     else
  3247.     {
  3248.       Textures[*Count] = Texture;
  3249.     }
  3250.  
  3251.     /* Test if this texture is already used. */
  3252.  
  3253.     for (i = 0; i < *Count; i++)
  3254.     {
  3255.       if (Textures[i] == Textures[*Count])
  3256.       {
  3257.         /* Add current weight to already existing texture weight. */
  3258.  
  3259.         Weights[i] += density;
  3260.  
  3261.         /* Any texture can only be in the list once --> exit. */
  3262.  
  3263.         return;
  3264.       }
  3265.     }
  3266.  
  3267.     Weights[(*Count)++] = density;
  3268.   }
  3269. }
  3270.  
  3271.  
  3272.  
  3273. /*****************************************************************************
  3274. *
  3275. * FUNCTION
  3276. *
  3277. *   Translate_Blob_Element
  3278. *
  3279. * INPUT
  3280. *
  3281. *   Element - Pointer to blob element
  3282. *   Vector  - Translation vector
  3283. *
  3284. * OUTPUT
  3285. *
  3286. *   Object
  3287. *
  3288. * RETURNS
  3289. *
  3290. * AUTHOR
  3291. *
  3292. *   Dieter Bayer
  3293. *
  3294. * DESCRIPTION
  3295. *
  3296. *   Translate a blob element.
  3297. *
  3298. * CHANGES
  3299. *
  3300. *   Sep 1994 : Creation.
  3301. *
  3302. ******************************************************************************/
  3303.  
  3304. void Translate_Blob_Element(Element, Vector)
  3305. BLOB_ELEMENT *Element;
  3306. VECTOR Vector;
  3307. {
  3308.   TRANSFORM Trans;
  3309.  
  3310.   if (Element->Trans == NULL)
  3311.   {
  3312.     /* This is a sphere component. */
  3313.  
  3314.     VAddEq(Element->O, Vector);
  3315.   }
  3316.   else
  3317.   {
  3318.     /* This is one of the other components. */
  3319.  
  3320.     Compute_Translation_Transform(&Trans, Vector);
  3321.  
  3322.     Transform_Blob_Element(Element, &Trans);
  3323.   }
  3324.  
  3325.   Translate_Textures(Element->Texture, &Trans);
  3326. }
  3327.  
  3328.  
  3329.  
  3330. /*****************************************************************************
  3331. *
  3332. * FUNCTION
  3333. *
  3334. *   Rotate_Blob_Element
  3335. *
  3336. * INPUT
  3337. *
  3338. *   Element - Pointer to blob element
  3339. *   Vector  - Translation vector
  3340. *
  3341. * OUTPUT
  3342. *
  3343. *   Object
  3344. *
  3345. * RETURNS
  3346. *
  3347. * AUTHOR
  3348. *
  3349. *   Dieter Bayer
  3350. *
  3351. * DESCRIPTION
  3352. *
  3353. *   Rotate a blob element.
  3354. *
  3355. * CHANGES
  3356. *
  3357. *   Sep 1994 : Creation.
  3358. *
  3359. ******************************************************************************/
  3360.  
  3361. void Rotate_Blob_Element(Element, Vector)
  3362. BLOB_ELEMENT *Element;
  3363. VECTOR Vector;
  3364. {
  3365.   TRANSFORM Trans;
  3366.  
  3367.   Compute_Rotation_Transform(&Trans, Vector);
  3368.  
  3369.   if (Element->Trans == NULL)
  3370.   {
  3371.     /* This is a sphere component. */
  3372.  
  3373.     MTransPoint(Element->O, Element->O, &Trans);
  3374.   }
  3375.   else
  3376.   {
  3377.     /* This is one of the other components. */
  3378.  
  3379.     Transform_Blob_Element(Element, &Trans);
  3380.   }
  3381.  
  3382.   Rotate_Textures(Element->Texture, &Trans);
  3383. }
  3384.  
  3385.  
  3386.  
  3387. /*****************************************************************************
  3388. *
  3389. * FUNCTION
  3390. *
  3391. *   Scale_Blob_Element
  3392. *
  3393. * INPUT
  3394. *
  3395. *   Element - Pointer to blob element
  3396. *   Vector  - Translation vector
  3397. *
  3398. * OUTPUT
  3399. *
  3400. *   Object
  3401. *
  3402. * RETURNS
  3403. *
  3404. * AUTHOR
  3405. *
  3406. *   Dieter Bayer
  3407. *
  3408. * DESCRIPTION
  3409. *
  3410. *   Scale a blob element.
  3411. *
  3412. * CHANGES
  3413. *
  3414. *   Sep 1994 : Creation.
  3415. *
  3416. ******************************************************************************/
  3417.  
  3418. void Scale_Blob_Element(Element, Vector)
  3419. BLOB_ELEMENT *Element;
  3420. VECTOR Vector;
  3421. {
  3422.   TRANSFORM Trans;
  3423.  
  3424.   if ((Vector[X] != Vector[Y]) || (Vector[X] != Vector[Z]))
  3425.   {
  3426.     if (Element->Trans == NULL)
  3427.     {
  3428.       /* This is a sphere component --> change to ellipsoid component. */
  3429.  
  3430.       Element->Type = BLOB_ELLIPSOID;
  3431.  
  3432.       Element->Trans = Create_Transform();
  3433.     }
  3434.   }
  3435.  
  3436.   if (Element->Trans == NULL)
  3437.   {
  3438.     /* This is a sphere component. */
  3439.  
  3440.     VScaleEq(Element->O, Vector[X]);
  3441.  
  3442.     Element->rad2 *= Sqr(Vector[X]);
  3443.   }
  3444.   else
  3445.   {
  3446.     /* This is one of the other components. */
  3447.  
  3448.     Compute_Scaling_Transform(&Trans, Vector);
  3449.  
  3450.     Transform_Blob_Element(Element, &Trans);
  3451.   }
  3452.  
  3453.   Scale_Textures(Element->Texture, &Trans);
  3454. }
  3455.  
  3456.  
  3457.  
  3458. /*****************************************************************************
  3459. *
  3460. * FUNCTION
  3461. *
  3462. *   Transform_Blob_Element
  3463. *
  3464. * INPUT
  3465. *
  3466. *   Element - Pointer to blob element
  3467. *   Trans   - Transformation
  3468. *
  3469. * OUTPUT
  3470. *
  3471. *   Object
  3472. *
  3473. * RETURNS
  3474. *
  3475. * AUTHOR
  3476. *
  3477. *   Dieter Bayer
  3478. *
  3479. * DESCRIPTION
  3480. *
  3481. *   Transform a blob element.
  3482. *
  3483. * CHANGES
  3484. *
  3485. *   Sep 1994 : Creation.
  3486. *
  3487. ******************************************************************************/
  3488.  
  3489. void Transform_Blob_Element(Element, Trans)
  3490. BLOB_ELEMENT *Element;
  3491. TRANSFORM *Trans;
  3492. {
  3493.   if (Element->Trans == NULL)
  3494.   {
  3495.     /* This is a sphere component --> change to ellipsoid component. */
  3496.  
  3497.     Element->Type = BLOB_ELLIPSOID;
  3498.  
  3499.     Element->Trans = Create_Transform();
  3500.   }
  3501.  
  3502.   Compose_Transforms(Element->Trans, Trans);
  3503.  
  3504.   Transform_Textures(Element->Texture, Trans);
  3505. }
  3506.  
  3507.  
  3508.  
  3509. /*****************************************************************************
  3510. *
  3511. * FUNCTION
  3512. *
  3513. *   Invert_Blob_Element
  3514. *
  3515. * INPUT
  3516. *
  3517. *   Element - Pointer to blob element
  3518. *
  3519. * OUTPUT
  3520. *
  3521. *   Object
  3522. *
  3523. * RETURNS
  3524. *
  3525. * AUTHOR
  3526. *
  3527. *   Dieter Bayer
  3528. *
  3529. * DESCRIPTION
  3530. *
  3531. *   Invert blob element by negating its strength.
  3532. *
  3533. * CHANGES
  3534. *
  3535. *   Sep 1994 : Creation.
  3536. *
  3537. ******************************************************************************/
  3538.  
  3539. void Invert_Blob_Element(Element)
  3540. BLOB_ELEMENT *Element;
  3541. {
  3542.   Element->c[2] *= -1.0;
  3543. }
  3544.  
  3545.  
  3546.  
  3547. /*****************************************************************************
  3548. *
  3549. * FUNCTION
  3550. *
  3551. *   Init_Blob_Queue
  3552. *
  3553. * INPUT
  3554. *
  3555. * OUTPUT
  3556. *   
  3557. * RETURNS
  3558. *   
  3559. * AUTHOR
  3560. *
  3561. *   Dieter Bayer
  3562. *   
  3563. * DESCRIPTION
  3564. *
  3565. *   Init queues for blob intersections.
  3566. *
  3567. * CHANGES
  3568. *
  3569. *   Dec 1994 : Creation.
  3570. *
  3571. ******************************************************************************/
  3572.  
  3573. void Init_Blob_Queue()
  3574. {
  3575.   Queue = (BSPHERE_TREE **)POV_MALLOC(Max_Queue_Size*sizeof(BSPHERE_TREE *), "blob queue");
  3576. }
  3577.  
  3578.  
  3579.  
  3580. /*****************************************************************************
  3581. *
  3582. * FUNCTION
  3583. *
  3584. *   Destroy_Blob_Queue
  3585. *
  3586. * INPUT
  3587. *
  3588. * OUTPUT
  3589. *
  3590. * RETURNS
  3591. *
  3592. * AUTHOR
  3593. *
  3594. *   Dieter Bayer
  3595. *
  3596. * DESCRIPTION
  3597. *
  3598. *   Destroy queues for blob intersections.
  3599. *
  3600. * CHANGES
  3601. *
  3602. *   Dec 1994 : Creation.
  3603. *
  3604. ******************************************************************************/
  3605.  
  3606. void Destroy_Blob_Queue()
  3607. {
  3608.   if (Queue != NULL)
  3609.   {
  3610.     POV_FREE(Queue);
  3611.   }
  3612.  
  3613.   Queue = NULL;
  3614. }
  3615.  
  3616.  
  3617.  
  3618. /*****************************************************************************
  3619. *
  3620. * FUNCTION
  3621. *
  3622. *   insert_node
  3623. *
  3624. * INPUT
  3625. *
  3626. * OUTPUT
  3627. *
  3628. * RETURNS
  3629. *
  3630. * AUTHOR
  3631. *
  3632. *   Dieter Bayer
  3633. *
  3634. * DESCRIPTION
  3635. *
  3636. *   Insert a node into the node queue.
  3637. *
  3638. * CHANGES
  3639. *
  3640. *   Feb 1995 : Creation.
  3641. *
  3642. ******************************************************************************/
  3643.  
  3644. static void insert_node(Node, size)
  3645. BSPHERE_TREE *Node;
  3646. int *size;
  3647. {
  3648.   /* Resize queue if necessary. */
  3649.  
  3650.   if (*size >= (int)Max_Queue_Size)
  3651.   {
  3652.     if (Max_Queue_Size >= INT_MAX/2)
  3653.     {
  3654.       Error("Blob queue overflow!\n");
  3655.     }
  3656.  
  3657.     Max_Queue_Size *= 2;
  3658.  
  3659.     Queue = (BSPHERE_TREE **)POV_REALLOC(Queue, Max_Queue_Size*sizeof(BSPHERE_TREE *), "blob queue");
  3660.   }
  3661.  
  3662.   Queue[(*size)++] = Node;
  3663. }
  3664.  
  3665.  
  3666.  
  3667. /*****************************************************************************
  3668. *
  3669. * FUNCTION
  3670. *
  3671. *   Create_Blob_Element_Texture_List
  3672. *
  3673. * INPUT
  3674. *
  3675. * OUTPUT
  3676. *
  3677. * RETURNS
  3678. *
  3679. * AUTHOR
  3680. *
  3681. *   Dieter Bayer
  3682. *
  3683. * DESCRIPTION
  3684. *
  3685. *   Create a list of all textures in the blob.
  3686. *
  3687. *   The list actually contains copies of the textures not
  3688. *   just references to them.
  3689. *
  3690. * CHANGES
  3691. *
  3692. *   Mar 1996 : Created.
  3693. *
  3694. ******************************************************************************/
  3695.  
  3696. void Create_Blob_Element_Texture_List(Blob, BlobList, npoints)
  3697. BLOB *Blob;
  3698. blobstackptr BlobList;
  3699. int npoints;
  3700. {
  3701.   int i, element_count, count;
  3702.   blobstackptr temp;
  3703.  
  3704.   if (npoints < 1)
  3705.   {
  3706.     Error("Need at least one component in a blob.");
  3707.   }
  3708.  
  3709.   /* Figure out how many components there will be. */
  3710.  
  3711.   temp = BlobList;
  3712.  
  3713.   for (i = count = 0; i < npoints; i++)
  3714.   {
  3715.     if (temp->elem.Type & BLOB_CYLINDER)
  3716.     {
  3717.       count += 3;
  3718.     }
  3719.     else
  3720.     {
  3721.       count++;
  3722.     }
  3723.  
  3724.     temp = temp->next;
  3725.  
  3726.     /* Test for too many components. [DB 12/94] */
  3727.  
  3728.     if (count >= MAX_BLOB_COMPONENTS)
  3729.     {
  3730.       Error("There are more than %d components in a blob.\n", MAX_BLOB_COMPONENTS);
  3731.     }
  3732.   }
  3733.  
  3734.   /* Allocate memory for components. */
  3735.  
  3736.   Blob->Data = (BLOB_DATA *)POV_MALLOC(sizeof(BLOB_DATA), "blob data");
  3737.  
  3738.   Blob->Data->References = 1;
  3739.  
  3740.   Blob->Data->Entry = (BLOB_ELEMENT *)POV_MALLOC(count*sizeof(BLOB_ELEMENT), "blob data");
  3741.  
  3742.   Blob->Data->Tree = NULL;
  3743.  
  3744.   /* Init components. */
  3745.  
  3746.   for (i = 0; i < count; i++)
  3747.   {
  3748.     init_blob_element(&Blob->Data->Entry[i]);
  3749.   }
  3750.  
  3751.   Blob->Data->Number_Of_Components = count;
  3752.  
  3753.   /* Allocate memory for list. */
  3754.  
  3755.   Blob->Element_Texture = (TEXTURE **)POV_MALLOC(count*sizeof(TEXTURE *), "blob texture list");
  3756.  
  3757.   for (i = 0; i < count; i++)
  3758.   {
  3759.     Blob->Element_Texture[i] = NULL;
  3760.   }
  3761.  
  3762.   for (i = element_count = 0; i < npoints; i++)
  3763.   {
  3764.     temp = BlobList;
  3765.  
  3766.     /* Copy textures. */
  3767.  
  3768.     switch (temp->elem.Type)
  3769.     {
  3770.       case BLOB_ELLIPSOID :
  3771.       case BLOB_SPHERE :
  3772.  
  3773.         /*
  3774.          * Copy texture into element texture list. This is neccessary
  3775.          * because individual textures have to be transformed too if
  3776.          * copies of the blob are transformed.
  3777.          */
  3778.  
  3779.         Blob->Element_Texture[element_count++] = Copy_Textures(temp->elem.Texture);
  3780.  
  3781.         break;
  3782.  
  3783.       case BLOB_CYLINDER :
  3784.  
  3785.         /*
  3786.          * Copy texture into element texture list. This is neccessary
  3787.          * because individual textures have to be transformed too if
  3788.          * copies of the blob are transformed.
  3789.          */
  3790.  
  3791.         Blob->Element_Texture[element_count++] = Copy_Textures(temp->elem.Texture);
  3792.  
  3793.         Blob->Element_Texture[element_count++] = Copy_Textures(temp->elem.Texture);
  3794.  
  3795.         Blob->Element_Texture[element_count++] = Copy_Textures(temp->elem.Texture);
  3796.  
  3797.         break;
  3798.     }
  3799.  
  3800.     BlobList = BlobList->next;
  3801.   }
  3802. }
  3803.  
  3804.  
  3805.  
  3806.